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Chapter  1 

Introduction  and  Literature  Review 


This  first  portion  of  this  chapter  describes  the  need  for  this  project  and  goals  upon  its  completion. 
The  second  portion  focuses  on  a  summary  of  previous  work  in  the  field. 

1.1  Objectives 

The  general  theory  for  the  removal  of  heat  from  concrete  by  embedded  cooling  pipes  was  first 
developed  for  use  in  the  design  of  the  Boulder  Canyon  (Hoover)  Dam  by  the  United  States 
Bureau  of  Reclamation  in  1949.  This  was  the  first  application  of  the  mathematical  model  of  heat 
conduction  into  a  format  that  could  be  utilized  in  the  design  and  implementation  of  cooling  coils 
to  cool  a  massive  concrete  system.  Subsequent  work  has  been  highly  dependent  on  the  theory 
developed  for  the  Hoover  Dam. 

1.1.1  Non-Linear  Incremental  Structural  Analysis 


Over  the  past  decade,  the  United  States  Army  Corp  of  Engineers  (USAGE)  has  developed  and 
implemented  the  use  of  nonlinear,  incremental  structural  analysis  (NISA)  procedures  to  predict 
the  effect  of  thermal  loads  due  to  the  heat  of  hydration  of  cement  in  massive  concrete  structures 
during  construction.  NISA’s  have  been  quite  successful  in  reducing  costs  of  construction  and/or 
providing  better  performance  through  the  construction  specifications.  When  a  dam  is 
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constructed,  it  is  built  in  sections  called  lifts  that  are  poured  one  at  a  time  in  a  particular  sequence 
determined  by  the  specifications  of  the  project.  NISAs  have  been  used  to  recommend  insulation 
requirements,  develop  appropriate  lift  dimensions,  specify  lift  placement  sequences,  and  specify 
concrete  constituents  in  order  to  reduce  construction  costs,  reduce  the  potential  for  cracking  and 
to  enhance  performance  of  the  structure. 

To  date,  the  commercial  program  ABAQUS  has  been  used  to  perform  a  majority  of  the  NISAs 
due  to  its  versatility.  The  ABAQUS  program  can  be  used  to  model  convection,  conduction,  heat 
generation,  cracking  criteria,  lift  placement,  age  dependent  material  properties  within  finite 
element  time  history  analyses.  Recently,  a  need  to  develop  a  modeling  procedure  for  NISAs  that 
can  account  for  cooling  coils  placed  within  massive  concrete  structures  has  arisen.  Often  the  heat 
generation  within  the  massive  concrete  structures  cannot  be  controlled  by  changing  the  concrete 
constituents,  reducing  lift  heights/widths,  or  modifying  the  construction  procedures.  Cooling 
coils  placed  within  the  concrete  are  needed  to  act  as  a  radiator,  constantly  carrying  heat  from  the 
source,  the  central  region  of  the  concrete,  in  order  to  reduce  the  thermal  gradient  within  the 
material. 

A  realistic  method  of  modeling  the  effects  of  cooling  coils  in  massive  concrete  structures  is 
presented  in  this  thesis  for  both  2D  and  3D  analyses.  Being  able  to  capture  the  thermal  changes 
in  the  concrete  and  cooling  coils  as  they  occur  in  time  is  the  primary  objective  of  the  modeling 
technique.  Previous  to  the  modeling  technique  presented  here,  there  was  no  acceptable  procedure 
for  modeling  cooling  coils  and  their  effects  within  massive  concrete  structures  using  ABAQUS. 

1.1.2  Portuguese  Dam  Project 

The  specific  motive  for  the  development  of  this  modeling  procedure  is  for  the  Portuguese  Dam. 
The  Portuguese  Dam  is  being  designed  by  the  Army  Corp  of  Engineers  at  the  time  of  this  study. 
The  modeling  procedure  developed  will  be  used  to  assist  in  the  design  of  the  cooling  coil  system. 
Parameters  from  the  Portuguese  Dam  project,  such  as  ambient  temperatures,  concrete  properties, 
and  water  properties  are  used  throughout  thesis  for  the  development  of  the  model.  However,  the 
model  can  be  used  for  any  massive  concrete  systems,  including  those  that  have  different 
parameters  and  specifications  than  the  Portuguese  Dam. 
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1.2  Organization 

Chapter  one  will  very  briefly  highlight  some  prior  work  involving  finite  element  analysis  of 
cooling  coils  embedded  in  mass  concrete  systems.  Chapter  two  will  describe  in  detail  the 
Boulder  Canyon  Dam  theory  for  the  removal  of  heat  by  embedded  pipes.  Chapter  two  will  also 
discuss  ABAQUS  capabilities  and  functions  and  the  preliminary  meshes  in  both  2D  and  3D  that 
were  used.  Chapter  3  will  discuss  the  requirements  of  the  model  and  modeling  procedures  that 
were  attempted.  It  analyzes  why  particular  methods  were  discarded.  Chapter  4  describes  in 
detail  the  modeling  procedures  used  to  model  the  embedded  cooling  coils,  which  is  the  use  of 
boundary  temperatures  nodes  in  ABAQUS.  Chapter  4  discusses  how  the  model  works, 
considerations,  concerns,  and  error.  Chapter  5  presents  results  obtained  from  using  the  modeling 
procedure.  The  results  presented  are  both  2D  and  3D,  and  include  2D  analysis  of  the  actual 
Portuguese  Dam  model.  Chapter  6  describes  how  to  implement  the  modeling  procedure  in  any 
mass  concrete  finite  element  model.  Chapter  7  contains  conclusions. 

1.3  Literature  Review 


Three  articles  involving  previous  work  using  finite  element  analysis  to  model  embedded  cooling 
coils  in  concrete  will  be  discussed  briefly  here.  In  Chapter  2,  the  Boulder  Canyon  Dam  theory  for 
embedded  cooling  coils  will  be  discussed  at  length. 

The  first  solution  to  the  problem  of  embedded  cooling  coils  used  to  remove  heat  from  mass 
concrete  systems  was  solved  by  the  Bureau  of  Reclamation  for  the  Boulder  Canyon  Dam  project 
in  1949.  The  theory  that  was  developed  had  certain  assumptions  that  will  be  discussed  further  in 
Chapter  2.  However,  it  produced  good  results  and  is  use  for  design  today.  Finite  element  models 
were  proposed  to  solve  the  problem,  such  as  that  by  E.L.  Wilson  at  the  University  of  Califomia- 
Berkeley  in  the  1960’s'.  These  models  were  extended  by  others  to  reduce  computational  time 


‘  Wilson,  Edward  L.  and  Robert  E.  Nickell.  ‘Application  of  the  Finite  Element  Method  to  Heat  Conduction 
Analysis.  ’  Nuclear  Engineering  and  Design.  North-Holland  Publishing  Company,  Amsterdam,  August 
1966. 
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using  a  substructure  methcxi  of  finite  element  analysis^  and  new  theory  has  been  developed^, 
including  that  for  nonmetal  pipes'*.  However,  until  this  point  there  was  not  a  finite  element 
modeling  procedure  available  using  ABAQUS  or  a  finite  element  modeling  procedure  that  could 
easily  be  used  by  designers  to  design  the  cooling  coils  systems  for  ntiassive  concrete  structures, 
particularly  for  NISA  analyses. 


^  Gong,  NG  and  GR  Li.  Stress  Analysis  in  Concrete  with  Cooling  Pipes.  Fifth  International  Conference  on 
Computing  in  Civil  and  Building  Engineering.  Anaheim,  California,  June  7-9,  1993. 

^  Liu,  X.-L.,  G.-R.  Li,  J.-S.  Jia,  and  P.  Xu.  A  Mathematical  Model  and  Finite  Element  Method  in  the 
Solution  of  Embedded  Pipe  Cooling  of  Concrete  Blocks.  Advanced  Computational  Methods  in  Heat 
Transfer  II,  Vol.l:  Conduction,  Radiation  and  Phase  Change.  Computational  Mechanics  Publications, 
Boston,  1992. 

“  Zhu,  Bofang.  ‘Effect  of  Cooling  by  Water  Flowing  in  NonMetal  Pipes  Embedded  in  Mass  Concrete.’ 
Journal  of  Construction  Engineering  and  Management.  Vol.  125,  No.  1,  January/February  1999. 
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Chapter  2 

Theory  and  Methodology 


2.1  Boulder  Canyon  Dam  Cement  and  Concrete 
Investigations  -  Cooling  of  Concrete  Dams 

Investigations  into  the  artificial  removal  of  heat  from  the  concrete  for  the  Hoover  Dam  project 
began  out  of  necessity.  When  massive  concrete  structures  are  constructed,  the  heat  of  hydration 
from  the  reaction  of  the  Portland  cement  and  water  in  the  concrete  is  released  and  thermal  stresses 
are  created  in  the  concrete.  These  thermal  stresses  can  cause  extensive  cracking  and  damage  to 
the  structural  system.  Due  to  the  massive  size  and  rapidity  of  construction  of  Hoover  Dam,  a 
major  design  issue  was  the  prevention  or  removal  of  this  excess  heat.  It  was  predicated  that, 
without  some  artificial  means  of  removing  the  excess  heat  of  hydration,  that  more  than  100  years 
would  be  required  for  the  concrete  in  the  dam  to  attain  thermal  and  volume  equilibrium^. 

2.1.1  Artificial  Cooling  by  Water  Flowing  in  Embedded  Pipes 


The  mathematical  theory  of  heat  conduction  was  used  to  develop  the  mathematical  model  for  the 
embedded  cooling  pipes  by  the  Bureau  of  Reclamation.  The  application  of  the  heat  conduction 
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theory  to  the  cooling  pipes  had  not  previously  been  documented  prior  to  the  construction  of 
Hoover  Dam.  The  Hoover  Dam  theory  has  since  been  used  by  designers  for  the  development  of 
cooling  coil  systems®.  The  basics  of  that  theory  will  be  outlined  here. 


To  begin  the  analysis,  consider  an  infinitely  long  hollow  circular  cylinder,  whose  inner  diameter 
is  ‘a’  and  outer  diameter  is  ‘b.’  The  cylinder  is  insulated  at  its  outer  boundary,  initially  at  a 
uniform  temperature,  and  has  zero  temperature  at  its  inner  surface.  The  governing  heat 
conduction  equation  is; 

dt  dr^  r  dr 

(2-1) 

The  variable  ‘r’  represents  the  radial  distance  to  any  point  in  the  cylinder  located  between  the 
outer  diameter  ‘b’  and  inner  diameter  ‘a.’  The  variable  6  represents  the  temperature,  h  is  the 
diffusivity,  and  t  is  time. 

The  mathematical  model  must  account  for  certain  boundary  conditions.  These  are  as  follows: 

1)  0=  0„  when  time  r  =  0 

2)  0  =  0  at  r  =  a  and  t-0 
un 

3)  —  =  0atr  =  bandt>0 
3r 

Assume  that  the  solution  is  of  the  form: 


6  =  dgUe 


(2-2) 


where  «  is  a  function  of  r  only,  a  is  a  constant,  and  h  is  the  diffusivity  of  the  concrete.  Using 
Bessel’s  differential  equations  of  order  zero,  a  solution  can  be  obtained  in  the  form  of  Eq.  2-2: 


®  United  States  Department  of  the  Interior,  Bureau  of  Reclamation.  Cooling  of  Concrete  Dams.  Boulder 
Canyon  Project.  Final  Reports.  Part  VII  -  Cement  and  Concrete  Investigations.  Bulletins,  1949.  Denver, 
CO,  page  2. 

®  For  a  more  thorough  explanation,  the  reader  should  see  United  States  Department  of  the  Interior,  Bureau 
of  Reclamation,  Cooling  of  Concrete  Dams.  Boulder  Canyon  Project,  Final  Reports.  Part  VII  -  Cement 
and  Concrete  Investigations.  Bulletin  3,  1949.  Denver,  CO. 
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-hWi 


/}=! 


(2-3) 


It  is  convenient  to  evaluate  Eq.  2-3  by  expressing  the  radius  a  as  some  constant  times  b  since  a 
definite  value  for  the  ratio  a/b  will  always  exist  for  any  given  hollow  cylinder.  If  b=D/2,  then 

—  =  y  (ab—)Ae 

(2-4) 

Therefore,  when  values  are  given  to  the  ratios  a/b  and  r/b,  the  ratio  & 6„  will  be  equal  for  two 
cylinders  of  different  sizes  and  thermal  properties  if  the  quantity  h^t/D^  has  the  same  value  for 
both  cylinders. 


To  obtain  the  mean  temperature  0,„  of  the  hollow  cylinder,  the  temperature  should  be  integrated 
throughout  the  volume  of  the  cylinder  and  divided  by  the  volume  of  the  cylinder  using  Eq.  2-3. 
The  resulting  expression  will  give  the  mean  temperature  at  a  time  t  for  an  infinite  hollow 
cylinder,  initially  at  a  uniform  temperature  6„  throughout,  whose  outer  boundary  is  insulated  and 
whose  inner  boundary  is  kept  at  zero  temperature.  After  introducing  a„b  and  b=D/2,  the 
following  expression  analogous  to  Eq.  2-4  results: 


00 


a. 


Ah‘{a„hyi 


lab  .  " 

-a'  h  a,b 


(2-5) 


By  inspection  of  Eq.  2-5,  it  is  apparent  that  if  a  definite  value  is  assigned  to  the  ratio  a/b,  the  ratio 
6J9o  will  have  the  same  value  for  two  cylinders  of  different  sizes  and  thermal  properties  if  the 
quantity  h^t/D^  is  equal  for  both  cylinders.  It  is  assumed  for  the  derivation  that  the  ratio  of  b/a  is 
equal  to  100.  Modifications  to  the  design  when  the  b/a  ratio  is  not  equal  to  100  will  be  discussed 
below. 


As  stated  above,  the  previous  discussion  was  completed  for  heat  flow  in  an  infinite  hollow 
cylinder,  initially  at  a  uniform  temperature  with  the  outer  boundary  perfectly  insulated  and  the 
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inner  boundary  held  at  zero  temperature  to  sirt^lify  the  calculations.  Therefore,  the  zero  point  of 
the  temperature  at  the  inner  boundary  must  be  shifted  to  the  cooling  coil  temperature.  The  results 
of  the  idealized  problem  of  the  infinite  hollow  cylinder  described  above  are  used  to  derive 
relations  for  the  actual  cooling  process  using  cooling  water  flowing  in  embedded  pipes. 

As  the  water  flows  through  a  cooling  pipe,  heat  is  transferred  from  the  concrete  to  the  water 
through  convection.  Therefore,  the  longer  the  water  is  traveling  through  the  cooling  coil  in  the 
concrete  block,  the  larger  its  temperature  increase.  This  results  in  different  rates  of  cooling  at 
different  locations  along  the  length  of  the  pipe  as  the  cooling  pipe  water  temperature  increases. 
The  concrete  closest  to  the  inlet  end  will  be  cooled  to  a  greater  extent  than  the  concrete  closer  to 
the  exit  of  the  pipe.  This  problem  has  been  addressed  by  the  Bureau  of  Reclamation  by 
specifying  that  the  cooling  coil  water  flow  be  reversed  every  twelve  hours  during  construction  of 
recent  dams,  leading  to  an  even  rate  of  cooling  in  the  concrete  for  design  purposes.  However,  the 
theory  created  for  Boulder  Dam  accounts  for  the  fact  that  different  cooling  rates  will  exist  for  all 
points  along  the  length  of  the  pipe. 


Solutions  to  the  actual  cooling  process  using  artificial  cooling  by  embedded  pipes  are  given  on 
three  sets  of  graphs  referred  to  as  the  X,  Y,  and  Z  curves.  All  three  graphs  depend  on  two 
dimensionless  parameters,  given  in  equations  2-6  and  2-7. 


KL 

h^t 


where: 

K  =  conductivity  of  concrete 
L  =  length  measure  along  the  cooling  pipe 
Cw  =  specific  heat  of  water  (or  other  cooling  fluid) 
pw  =  density  of  water  (or  other  cooling  fluid) 

Q  =  cooling  fluid  flow  rate 
h  =  diffusivity  of  concrete 


(2-6) 


(2-7) 


t  =  time  since  cooling  commenced 
D  =  diameter  of  the  cooled  cylinder 


9 


The  units  must  be  kept  consistent  in  order  to  create  dimensionless  parameters.  The  three  graphs 
are  defined  as  follows: 


X  = 


Mean  temperature  of  concrete  cylinder  of  length  L  —  inital  temperature  of  water 
Initial  temperature  of  concrete  —  initial  temperature  of  water 


Temperature  of  water  at  given  distance  -  inital  temperature  of  water 
Initial  temperature  of  concrete  -  initial  temperature  of  water 


Mean  temperature  of  concrete  at  a  given  length  L  -  inital  temperature  of  water 
Initial  temperature  of  concrete  -  initial  temperature  of  water 


The  Y-curves  are  used  to  determine  the  temperature  rise  of  the  cooling  coil  water  as  it  travels 
through  the  concrete.  The  X-curves  represent  the  final  mean  temperature  difference  in  degrees 
per  degree  of  initial  temperature  difference.  The  Z-curves  are  used  to  compute  the  mean 
temperature  of  concrete  at  a  given  length  from  the  inlet.  As  stated  earlier,  it  is  current  procedure 
of  the  Bureau  of  Reclamation  to  require  that  the  flow  of  cooling  coil  water  be  reversed  every 
twelve  hours.  As  a  result,  the  mean  temperature  of  the  concrete  is  assumed  to  be  relatively 
constant  throughout  the  block  and  the  Z-curves  are  rarely  used. 


2.1.2  Development  of  Y-Curves 


Across  the  surface  r=a,  the  heat  H,,  which  flows  out  of  the  infinite  hollow  cylinder  per  unit  length 
of  cylinder  per  unit  time  is: 


=  Tcif^  -a^)cp 


de 

m 

dt 


(2-8) 


By  differentiating  Eq.  2-5  with  respect  to  time  (r),  we  obtain  the  rate  of  change  of  the  mean 
cylinder  temperature  0m  with  respect  to  time: 
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Sab 


h^0„ 


u2  „2 

b  —a 


Sulfanb7V„(a„b> 


n=l 


4hMgnbyi 


(2-9) 


Substituting  k-h^cp  and  Eq.  2-9  into  2-8,  we  obtain: 


H..  =  kOR 


where  R  is  defined  in  Eq.  2-11: 


R 


4a7C^  Y  , 


D 


4h^(an>>y» 


E^olM-g- K(a„b> 


(2-10) 


(2-11) 


As  a  result  of  Eq.  2-11,  after  defining  a/b  and  the  roots  ( a„b),  values  of  R  may  be  determined  as  a 
function  of  h^t/D^.  When  h^t/D^=0,  R  is  equal  to  R„. 


As  stated  earlier,  the  development  of  this  theory  accounts  for  the  fact  that  the  water  temperature 
rises  as  it  flows  through  the  concrete  cylinder  and  this  accounts  for  varying  rates  of  heat  flow 
from  the  concrete  to  the  water  in  the  actual  cylinder.  In  order  to  do  this,  a  concrete  cylinder  of 
finite  length  is  considered.  At  a  distance  L  from  the  inlet  end  of  the  cylinder  and  a  time  after 
cooling  has  begun,  a  temperature  differential  in  the  cooling  water  of  6„(dY/d^)dX  will  occur. 
Substituting  this  temperature  differential  into  Eq.  2-10  for  by  making  the  appropriate  change  in 
R(h^t/D^)  and  integrating  from  0  to  r  yields  a  heat  flow  per  unit  length  of  cylinder  per  unit  time. 
Summing  up  these  temperature  differentials  over  the  length  of  the  cylinder  by  performing  a 
second  integration  with  respect  to  L  gives  the  total  heat  flow  per  unit  of  time: 


(t-x) 


hn 


D' 


^d?idL 

dX 


(2-12) 
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Eq.  2-12  is  for  a  hollow  concrete  cylinder  of  length  L,  with  initial  conditions  as  stated  earlier. 

This  heat  will  be  absorbed  by  the  water  in  the  cooling  pipe.  The  amount  of  heat  absorbed  per  unit 
time  is: 


(2-13) 

where  the  w  subscript  signifies  that  the  variable  is  a  property  of  water  and  signifies  the 
temperature  rise  of  the  water  at  a  given  time  and  distance  from  the  inlet  end.  If  the  incoming 
water  temperature  is  equal  to  0,  then  Y=6je„.  Substituting  this  expression  in  Eq.  2-13  gives: 


^7  =  c^p„q„e„Y 


(2-14) 


Equating  equations  2-14  and  2-12  gives: 


c.P.q.6.Y  =  Ke,0R 


(.-X) 


hh 


D' 


(2-15) 


Eq.  2-15  is  somewhat  difficult  to  solve  because  it  includes  Y  under  an  integral  on  one  side  and 
independently  on  the  other  side.  For  this  reason,  the  equation  will  first  be  solved  for  the  special 
case  of  h^t/D^=0.  In  this  situation,  the  temperature  of  the  concrete  is  d,„  the  temperature  of  the 
water  at  any  point  along  the  cylinder  is  and  R  is  equal  to  /?„as  discussed  earlier.  Under  these 
conditions,  Eq.  2-15  becomes: 

L 

c.P.qAY  =  Ke.Rj(l-Y)dL 

0 

(2-16) 


The  time  the  cooling  water  takes  to  flow  through  any  given  coil  is  assumed  to  be  small  compared 
with  the  time  necessary  to  cool  the  surrounding  concrete.  Let: 


KL 


^wPwqw 


= - dL 

^wPw^w 


(2-17) 
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Substitute  Eq.  2-17  into  Eq.  2-16  and  obtain: 

KL 

Y  =  R„  ” j(l-Y)d^ 

0 

(2-18) 


If  Eq.  2-18  is  differentiated  with  respect  to  we  obtain: 

dY 


or 


dY 


(2-19) 


This  is  a  first  order  linear  differential  equation  which  has  a  solution  of  the  form: 

(2-20) 


When  ^  and  Y  are  equal  to  0,  and  constant  C/  is  equal  to  -1,  then  Eq.  2-20  becomes: 


Y  =  l-e"''“^ 


(2-21) 


The  value  of  Ro  is  determined  from  Eq.  2-11  with  h^t/D^  equal  to  0  and  assuming  that  the  ratio  of 
b/a=l(y0.  If  the  first  five  values  of  the  series  in  Eq.  2-11  are  used,  then  Ro  is  equal  to  2.693590. 
Substituting  this  value  in  Eq.  2-21  obtains: 


Y  =  l-e 


(2-22) 


Eq.  2-22  is  used  to  plot  the  curve  for  h^t/D^  equal  to  0  as  shown  in  Figure  2-1. 
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Solving  Eq.  2-14  for  the  case  of  /zVD^  0  is  a  more  complicated  process  since  Y  is  included 
within  the  integral  and  also  independently  on  the  other  side  of  the  equation  as  stated  earlier.  To 
begin,  let: 


o=|r 


(t-^) 


D' 


dX 


(2-22) 


Substituting  Eq.  2-22  into  Eq.  2-14  obtains: 


CwPw 


(2-23) 


Using  the  substitutions  given  in  Eq.  1-17,  Eq.  1-21  becomes: 

KL 

^wPwQw 

Y=  jGd^ 

(2-24) 

Since  the  above  equations  include  Y  independently  and  under  an  integral,  the  solution  must 
follow  an  iterative  procedure  that  will  not  be  described  in  detail  here^.  Using  this  iterative 
procedure,  the  Y-curves  for  values  of  h^t/D^  greater  than  0  are  obtained.  The  Y-curves  are  shown 
in  Figure  2-1. 


’  For  a  description  of  this  iterative  procedure,  please  see  United  States  Department  of  the  Interior,  Bureau 
of  Reclamation.  Cooling  of  Concrete  Dams.  Boulder  Canyon  Project.  Final  Reports.  Part  VII  -  Cement 
and  Concrete  Investigations.  Bulletins,  1949.  Denver,  CO.  Pages  126-127. 
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>- 


Figure  2-1;  Curves  ofYvs.  KUc„pji„for  Various  Values  of  hH/D^^ 


*  United  States  Department  of  the  Interior,  Bureau  of  Reclamation.  Cooling  of  Concrete  Dams.  Boulder 
Canyon  Project.  Final  Reports.  Part  VII  -  Cement  and  Concrete  Investigations.  Bulletin  3,  1949.  Denver, 


CO. 


2.1.3  Development  of  Z-Curves 
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The  Z-curves  allow  us  to  obtain  the  mean  temperature  of  the  concrete  cylinder  at  any  section 
along  the  cooling  pipe.  This  mean  concrete  temperature  can  be  obtained  by  summing  all  the 
increments  of  mean  concrete  temperature  produced  by  changes  in  the  temperature  of  the  cooling 
water  up  to  the  time  t,  added  to  the  temperature  of  the  cooling  water  at  time  t.  An  increment  of 
mean  temperature  which  occurs  at  time  X  will  have  been  in  existence  for  a  period  at  time  t. 
The  increment  of  mean  temperature  produced  by  an  incremental  change  of  6„(dY/dX)dX  in  the 
temperature  of  the  cooling  water  may  be  obtained  by  Eq.  2-25: 


4h^(a„b)Tl->-) 


(dej. 


2ab0,  ay 
b'-a'  dX 


d^X- 


u-(a„b-)A„e 


n=l 


ttnb 


(2-25) 


In  functional  notation,  Eq.  2-25  can  be  expressed: 


where: 


(deJi=e„- 


00 


(t-^> 


ill 

D' 


ay 


dx 


00 


2ab 


u. 


a  b- 


A  e 


a„b 


(2-26) 


The  sum  of  all  these  increments  up  to  the  time  t  will  be: 


(2-27) 
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Then,  if  Z  is  the  ratio  of  the  mean  temperature  at  the  time  t  of  any  section  along  the  length  of  the 
cylinder  to  the  temperature  6„,  then: 

e.z=e,Y+(eJ, 

(2-28) 


Combining  Eq.  2-25  and  Eq.  2-26: 

Z=Y-t-J 


0. 


(t-^> 


D' 


dX 


dX 


(2-29) 


Eq.  2-29  is  again  solved  by  graphical  iteration  that  will  not  be  discussed  here®.  This  process 
obtains  the  Z-curves  as  shown  in  Figure  2-2. 


2.1.4  Development  of  X-Curves 


The  X-curves  are  used  to  obtain  the  mean  temperature  of  the  concrete  cylinder  along  its  entire 
length  from  the  inlet  end  to  any  particular  section  along  the  length.  This  value  can  be  obtained  by 
taking  the  Z-curve  from  zero  to  the  length  chosen  and  dividing  by  the  length.  This  is  described 
mathematically  by: 


(2-30) 


Using  the  variable  substitutions  from  Eq.  2-17,  the  X-curves  are  obtained  as  shown  in  Figure  2-3. 

KL 


_  ^wPy/lw 

KL 


Jzaf 


(2-31) 


’  For  a  description  of  this  iterative  procedure,  please  see  United  States  Department  of  the  Interior,  Bureau 
of  Reclamation.  Cooling  of  Concrete  Dams.  Boulder  Canyon  Project.  Final  Reports.  Part  VII  -  Cement 
and  Concrete  Investigations.  Bulletin  3, 1949.  Denver,  CO.  Page  129. 
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Figure  2-2:  Curves  ofZ  vs.  KUc^pji^for  Various  Values  ofhH/D^‘° 


*°  United  States  Department  of  the  Interior,  Bureau  of  Reclamation.  Cooling  of  Concrete  Dams.  Boulder 
Canyon  Project.  Final  Reports.  Part  VII  -  Cement  and  Concrete  Investigations.  Bulletins,  1949.  Denver, 
CO. 


Figure  2-3;  Curves  ofX  vi'.  KL/c^p^„for  Various  Values  ofh^t/D^" 


”  United  States  Department  of  the  Interior,  Bureau  of  Reclamation.  Cooling  of  Concrete  Dams.  Boulder 
Canyon  Project.  Final  Reports.  Part  VII  -  Cement  and  Concrete  Investigations.  Bulletin  3,  1949.  Denver, 
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It  is  important  to  recognize  at  this  point  the  development  of  the  X,  Y,  and  Z  curves  assumed  a  b/a 
ratio  equal  to  100,  as  stated  earlier.  However,  a  b/a  ratio  equal  to  100  rarely  occurs  in  practice. 

To  account  for  this  discrepancy,  but  still  allow  for  use  of  the  X,  Y,  and  Z  curves,  a  modification  is 
made  to  both  the  diffusivity  value  (/i)  and  the  diameter  of  the  cooled  cylinder  (D).  These 
modifications  have  been  calculated  by  the  Bureau  of  Reclamation  for  common  lift  dimensions,  as 
shown  in  Table  2-1. 

Table  2-1:  Values  ofD,  D^,  and  h^/for  Pipe  Cooling‘s 


D 

Vertical 

Horizontal 

n  f 

feet 

feet 

feet 

1^31 

ft"/hr 

2  1/2 

2  1/2 

2.82 

7.95 

1.31  xh^ 

5 

2  1/2 

3.99 

in  kmH 

1.19 

5 

3 

4.35 

1.16 

5 

4 

5.02 

25.2 

1.12 

5 

5 

5.64 

31.81 

1.09 

5 

6 

6.18 

38.19 

1.07 

7  1/2 

2  1/2 

4.88 

23.81 

1.13  xh^ 

7  1/2 

4 

6.15 

37.82 

1.07 

7  1/2 

5 

6.86 

47.06 

1.04 

7  1/2 

6 

7.54 

56.85 

1.02 

7  1/2 

7  1/2 

8.46 

71.57 

1 

7  1/2 

9 

9.26 

85.75 

0.98 

10 

2  1/2 

5.64 

31.83 

1.09  xh" 

10 

4 

7.14 

50.93 

1.03 

10 

5 

7.98 

.  63.66 

1.01 

10 

6 

8.74 

76.39 

0.99 

10 

7 

9.44 

89.13 

0.97 

10 

8 

10.09 

101.86 

0.96 

10 

■■ 

10.7 

114.59 

0.95 

10 

11.28 

127.32 

0.94 

10 

12.36 

152.79 

0.92 

10 

13.35 

178.25 

0.91 

10 

13.82 

190.99 

0.9 

From  Control  ofCracldng  in  Mass  Concrete  Structures.  C.L.  Townsend,  United  States  Department  of 
the  Interior,  Bureau  of  Reclamation.  October  1965.  Denver,  CO. 
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Another  important  point  worth  noting  to  is  that  the  time  value  used  for  the  h^t/lf  term  on  the 
curves  is  the  total  time  after  the  start  of  cooling  and  includes  no  heat  of  hydration  of  the  concrete 
system.  The  consequences  of  this  will  be  discussed  extensively  in  Section  4.1.3. 


2.2  ABAQUS 


ABAQUS  is  an  engineering  modeling  program  that  solves  problems  using  finite  element  analysis. 
ABAQUS  can  solve  a  wide  range  of  problems,  from  simple  linear  analyses  to  complex  non-linear 
problems,  including  stress/displacement,  heat  transfer,  mass  diffusion,  thermal  management  of 
electrical  components,  acoustics,  soil  mechanics,  and  piezoelectric  analyses.  ABAQUS  has 
extensive  options  available  for  the  element  shapes,  dimensions,  and  material  types.  The  modeling 
procedure  discussed  in  this  paper  uses  ABAQUS/standard,  which  is  a  general-purpose  analysis 
module  that  can  solve  linear  and  nonlinear  static,  dynamic,  thermal,  and  electrical  problems. 

ABAQUS  also  offers  the  incorporation  of  user  written  subroutines.  Subroutines  are  offered  as 
options  for  user  control  of  variables  for  particular  problems.  For  example,  if  heat  generation  is 
required  for  the  material  type  modeled  in  the  analysis,  ABAQUS  can  incorporate  a  user 
subroutine  that  contains  the  FORTRAN  code  for  the  calculation  of  the  heat  generated  during  the 
analysis.  ABAQUS  then  incorporates  this  subroutine  into  its  analysis.  User  subroutines  can  only 
be  used  if  ABAQUS  provides  the  option  of  having  user  control  of  a  variable.  They  are  provided 
to  increase  the  functionality  of  the  ABAQUS  options  for  which  data  line  usage  or  predetermined 
variable  values  alone  may  be  too  restrictive. 

2.2.1  Heat  Transfer  Analysis 


ABAQUS  can  solve  uncoupled  heat  transfer  analyses,  sequentially  coupled  thermal-stress 
analyses,  fully  coupled  thermal-stress  analyses,  adiabatic  analyses,  coupled  thermal-electrical, 
and  cavity  radiation  problems.  For  the  modeling  of  the  cooling  coils,  an  uncoupled  heat  transfer 
analysis  was  used.  After  obtaining  the  temperature  distribution  from  the  uncoupled  heat  transfer 
analysis,  these  temperatures  can  be  used  to  analyze  the  thermal  stress  in  the  system. 
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Uncoupled  heat  transfer  analysis  is  used  to  model  solid  body  heat  conduction  with  general, 
temperature  dependent  conductivity,  internal  energy,  and  general  convection  and  radiation 
conditions.  Forced  convection  of  a  fluid  through  a  mesh  can  be  modeled  using  specified 
elements.  This  will  be  discussed  further  in  section  3.5.  These  problems  can  be  transient  or 
steady-state,  linear  or  nonlinear,  and  can  include  thermal  interactions  such  as  gap  radiation, 
conductance,  and  heat  generation. 

Completing  a  transient  heat  transfer  analysis  with  second  order  elements  in  ABAQUS  can  lead  to 
the  appearance  of  spurious  oscillations  in  the  solution,  particularly  in  the  vicinity  of  boundaries 
with  rapidly  changing  temperatures.  To  avoid  this  inaccuracy,  the  characteristic  length  of  the 
element  and  the  time  increment  should  satisfy  the  relationship: 

(2-32) 

6k 

where  At  is  the  time  increment,  p  is  the  density,  c  is  the  specific  heat ,  k  is  the  thermal 
conductivity,  and  Al  is  the  characteristic  length  of  the  element  (in  this  case  the  distance  between 
two  nodes). 

2.3  Development  of  a  Preliminary  Model 

The  development  of  this  model  was  based  on  the  parameters  used  to  construct  the  Portuguese 
Dam  in  Puerto  Rico.  The  Army  Corp  of  Engineers  specifications  state  that  the  placement 
temperature  of  the  concrete  for  this  project  is  75°F.  Heat  of  hydration  of  the  Portland  Cement 
reaction  with  water  is  calculated  and  input  in  the  ABAQUS  model  using  the  subroutine 
HETVAL.  At  the  end  of  each  increment,  the  HETVAL  subroutine  is  called.  This  routine 
calculates  the  additional  heat  added  to  the  concrete  due  to  the  heat  of  hydration  based  on  the 
length  of  time  the  particular  lift  has  been  in  place.  This  additional  heat  is  then  accounted  for  by 
ABAQUS.  HETVAL  will  calculate  different  values  depending  on  the  lift  location  and  time  since 
its  placement.  After  28  days,  the  heat  of  hydration  is  assumed  to  be  negligible. 

ABAQUS  allows  for  the  input  of  boundary  conditions  to  be  specified  for  the  analysis.  These 
boundary  conditions  are  prescribed  at  nodes  for  specified  degrees  of  freedom,  such  as 
displacement  or  temperature.  These  boundary  conditions  can  be  constant,  or  they  can  be  defined 
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to  vary  according  to  an  amplitude  curve  that  defines  the  magnitude  of  the  prescribed  boundary 
conditions.  The  lifts  will  be  cooled  (or  heated)  by  external  temperature  conditions  as  the  exterior 
faces  of  the  lifts  are  exposed  to  the  ambient  temperature  of  the  surrounding  air.  For  this  particular 
project,  ambient  temperature  conditions  at  the  dam  site  in  Puerto  Rico  were  used  as  boundary 
temperatures.  The  ambient  temperature  was  defined  by  an  amplitude  function  defining  the 
seasonal  variation  in  temperature,  with  an  average  temperature  of  78°F,  with  83°F  and  73°F  as  the 
two  extremes. 

Knowledge  of  the  cooling  pipe  characteristics  was  necessary  to  develop  the  model.  The  cooling 
fluid  used  in  the  model  is  water  with  a  flow  rate  of  four  gallons  per  minute.  The  cooling  coils 
themselves  had  a  one  inch  outer  diameter  and  a  0.887  inch  inner  diameter.  The  density  of 
concrete  is  156.3  psi,  conductivity  is  19.131 12  BTU/day/°F/ft.,  and  specific  heat  is  0.17 
BTU/lb/°F.  It  is  important  to  note  that  none  of  these  values  are  set;  the  model  allows  for  the 
variability  of  any  of  these  parameters. 

2.3.1  2D  Preliminary  Mesh 

A  2D  preliminary  mesh  was  developed  to  test  various  methods  of  modeling  the  cooling  coils. 

This  mesh  was  small  and  simple  to  ensure  that  an  attempted  method  could  be  analyzed  without 
complications  from  the  attributes  of  the  mesh.  The  Army  Corp  of  Engineers  uses  2.5-foot 
element  dimensions  for  the  Portuguese  Dam,  determined  by  the  material  properties  as  shown  in 
Eq.  2-32.  The  preliminary  mesh  developed  here  used  two-foot  elements  for  two  reasons.  The 
first  was  to  ensure  that  Eq.  2-32  was  satisfied.  Secondly,  it  was  anticipated  that  the  cooling  coils 
may  need  to  be  placed  at  node  locations.  If  the  total  length  of  the  elements  was  two  feet,  this 
allowed  for  the  cooling  coils  to  be  placed  at  various  spacings  that  are  multiples  of  two  feet. 

When  a  dam  is  constructed,  it  is  built  in  sections  called  lifts  that  are  poured  one  at  a  time  in  a 
particular  sequence  determined  by  the  specifications  of  the  project.  The  lift  dimensions  are 
typically  five  to  ten  feet  in  the  vertical  direction,  with  varying  horizontal  dimensions.  For  the 
preliminary  mesh,  two  lifts  were  modeled.  Each  lift  had  a  vertical  height  of  ten  feet  and  a 
horizontal  dimension  of  16  feet.  During  construction,  vertical  lifts  will  be  placed  in  seven-day 
increments.  As  a  lift  is  finished,  cooling  colls  are  laid  on  the  top  of  the  finished  lift  and  the  new 
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lift  will  be  poured  directly  on  top  of  the  cooling  coils.  This  sequencing  means  that  each  lift  will 
be  cooled  by  cooling  coils  at  the  boundaries  between  the  lift  interface  above  and  below  the 
current  lift.  Therefore,  the  two  lift  preliminary  mesh  was  analyzed  for  cooling  coils  placed  at  the 
two-lift  interface. 

Elements  used  in  the  model  were  second  order  heat  transfer  elements  type  DC2D8.  These  are 
eight  node  quadratic  elements  with  nine  integration  points  as  shown  in  Figure  2-4. 

Figure  2-4:  8-Node  Quadratic  Element 
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Second-order  elements  provide  higher  accuracy  than  first-order  elements  for  smooth  problems 
that  do  not  involve  complex  contact  conditions,  impact,  or  severe  element  distortions  and  reduce 
the  number  of  elements  required.  The  DC2D8  elements  have  an  active  degree  of  freedom  ‘11’ 
which  is  the  temperature  degree  of  freedom.  Material  characteristics  were  defined  for  all 
elements  with  the  concrete  parameters  discussed  above.  The  2D  preliminary  mesh  is  shown  in 
Figure  2-5.  The  ambient  temperature  was  applied  as  the  top  and  bottom  boundary  condition,  and 
symmetry  boundaries  were  applied  to  the  right  and  left  side. 


Figure  2-5:  2D  Preliminary  Mesh 
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2.3.2  3D  Preliminary  Mesh 


As  in  the  2D  mesh,  the  3D  mesh  included  two,  ten-foot  vertical  lifts.  The  lifts  had  a  horizontal 
width  of  16  feet,  and  a  depth  of  10  feet.  Again  the  cooling  coils  were  placed  at  the  two-lift 
interface.  Element  type  DC3D20  was  used.  These  elements  are  20-node  quadratic  brick 
elements  with  an  active  degree  of  freedom  ‘11’  which  is  the  temperature  degree  of  freedom.  This 
element  is  shown  in  Figure  2-6. 

Figure  2-6:  20-node  Element 
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Figure  2-7  shows  the  integration  scheme  in  the  layer  closest  to  the  1 -2-3-4  faces.  The  integration 


points  in  the  second  and  third  layers  are  numbered  consecutively. 
Figure  2-7:  20-node  Element  Integration  Points 


The  3D  preliminary  mesh  is  shown  in  Figure  2-8. 
Figure  2-8;  3D  Preliminary  Mesh 


2.4  Summary 
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The  theory  developed  to  design  the  cooling  system  for  Boulder  Canyon  Dam  by  the  Bureau  of 
Reclamation  has  been  used  in  the  past  with  good  results  to  predict  the  concrete  temperatures  in  a 
massive  concrete  structure  cooled  by  embedded  coils.  This  theory  will  prove  valuable  to  the 
development  of  the  model  presented  in  this  thesis  as  well.  The  preliminary  meshes,  both  in  2D 
and  3D,  presented  in  this  chapter  will  be  used  in  Chapter  3  to  evaluate  the  validity  of  potential 
modeling  procedures. 
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Chapter  3 

Modeling  Procedures 


This  chapter  will  outline  the  general  approach  taken  to  determine  the  best  method  for  modeling 
the  cooling  coils.  Approaches  were  attempted  and  evaluated  for  validity  using  the  simple 
preliminary  meshes  described  in  Chapter  2.  Using  the  preliminary  mesh,  the  effect  of  the 
methods  attempted  could  be  easily  evaluated  because  the  mesh  itself  was  uncomplicated. 

3.1  Model  Requirements 

The  Army  Corp  of  Engineers  wanted  to  implement  a  procedure  to  be  used  to  evaluate  cooling  coil 
designs  in  both  two  dimensional  and  three-dimensional  models.  Anticipating  the  time  of  analysis 
and  complexity  required  for  the  3D  model,  the  2D  model  will  be  used  to  evaluate  different 
cooling  coil  spacings.  When  the  best  spacing  is  determined,  this  spacing  will  be  implemented  in 
a  3D  model  to  find  the  temperature  distribution  in  the  concrete.  These  temperatures  will  then  be 
used  to  determine  the  thermal  stresses.  The  procedures  used  to  evaluate  the  cooling  coils  could 
be  different  for  both  the  2D  and  3D.  However,  if  they  were  the  same  or  similar,  it  would  help 
minimize  confusion  and  mistakes  during  implementation. 

The  procedure  used  to  implement  the  cooling  coils  in  an  ABAQUS  model  had  several 
requirements.  The  procedure  itself  was  to  be  used  by  engineers  at  the  Army  Corp  of  Engineers 
for  the  design  phase  of  massive  concrete  systems  to  help  design  the  cooling  coils.  Typically,  the 
design  engineers  will  have  already  created  complicated  2D  and  3D  models  of  the  proposed  dam 
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or  other  structure.  Therefore,  the  cooling  coil  procedure  needed  to  be  easily  input  into  an  existing 
model  without  modifying  the  major  characteristics  of  the  existing  model  such  as  the  element  size 
or  node  spacing.  In  addition,  the  procedure  needed  to  be  as  uncomplicated  as  possible  to 
minimize  the  time  required  for  implementation  into  an  existing  model  and  to  limit  the  frequency 
of  errors  made. 

It  was  also  necessary  to  limit  modifications  to  the  procedure  that  are  dependent  on  the  specific 
problem  being  designed.  The  more  changes  the  user  would  need  to  make  to  the  procedure  or 
code,  the  more  opportunities  for  errors  during  use,  particularly  if  these  changes  required 
FORTRAN  coding  by  the  user. 

The  amount  of  heat  that  the  cooling  coil  removes  from  the  concrete  is  dependent  on  the 
differential  temperature  between  the  concrete  and  the  water.  The  input  temperature  of  the  water 
is  constant  because  the  construction  specifications  require  that  the  water  leaving  the  concrete  be 
cooled  to  the  specified  inlet  temperature  before  being  recycled  through  the  concrete.  As  the 
cooling  water  moves  through  the  concrete,  it  becomes  warmer  as  it  absorbs  heat  from  the 
concrete.  Therefore,  the  concrete  at  the  input  point  of  the  coil  will  have  a  larger  rate  of  cooling 
than  the  concrete  at  the  exit  point  of  the  water  because  the  temperature  differential  will  be  larger 
at  the  input.  To  avoid  the  occurrence  of  different  rates  of  cooling  in  different  portions  of  the  lift, 
it  is  assumed  that  construction  specifications  require  that  the  water  flow  be  reversed  every  twelve 
hours.  Therefore  it  will  be  assumed  that  the  rate  of  cooling  through  the  concrete  block  is 
relatively  constant. 

The  last  concern  is  that  the  cooling  coils  are  not  operational  at  all  times.  They  are  used  initially  to 
control  the  peak  temperature  due  to  heat  generation.  Then  they  are  shut  off,  and  turned  on  again 
to  help  control  the  cooling  regulating  the  closure  of  the  grouted  joints.  Therefore,  a  method  that 
could  be  turned  off  and  on  at  different  points  in  the  analysis  was  necessary. 

3.2  CFLUX  and  DFLUX 

The  first  method  attempted  was  using  the  two  flux  options  available  from  ABAQUS.  The  flux 
options  are  a  type  of  thermal  loading.  CFLUX  is  used  to  specify  a  concentrated  heat  flux  at  a 
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node.  The  DFLUX  option  is  used  to  specify  distributed  surface  fluxes  (on  element  faces)  or  body 
fluxes  (flux  per  unit  volume).  The  flux  option  can  be  used  to  specify  the  rate  at  which  heat  is 
leaving  the  mesh.  Using  DFLUX  is  not  an  option  because  the  cooling  coils  are  very  small  (1” 
diameter)  compared  to  the  element  surface  (2’  x  2’).  Therefore,  applying  a  heat  flux  over  the 
element  surface  would  not  be  an  accurate  representation  of  the  manner  in  which  the  cooling  coil 
removes  heat  from  the  concrete. 

Using  CFLUX,  an  arbitrary  heat  flux  was  prescribed  at  a  node  to  determine  whether  this  would 
remove  heat  from  the  concrete  mesh  as  shown  in  Figure  3-1.  All  nodes  had  an  initial  temperature 
of  100°F  and  heat  generation  was  included.  The  length  of  the  run  was  2.5  days. 

Figure  3-1:  CFLUX  (Scale:  maximum  145°F,  minimum  92. FF,  interval  4‘’F) 


The  cooled  edges  on  the  top  and  bottom  of  the  mesh  are  the  boundaries  with  the  ambient 
temperature  applied.  This  plot  demonstrates  that  applying  a  flux  to  the  central  node  in  the  mesh 
did  remove  heat  from  the  system  via  the  central  node  where  the  flux  was  applied.  However, 
CFLUX  does  not  have  a  user  subroutine  available  to  allow  the  user  to  write  code  to  modify  the 
flux  rate.  It  only  allows  for  the  flux  rate  to  be  constant  or  controlled  by  a  prescribed  amplitude 
function.  This  is  not  adequate  because  the  rate  of  cooling  will  change  with  time  since  due  to  the 
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effects  of  heat  generation  and  the  difference  in  the  heat  absorbed  by  the  water  as  the  temperature 
of  the  concrete  increases  and  decreases. 

3.3  Gap  Elements 

Gap  elements  allow  for  contact  between  two  nodes  which  are  either  adjacent  to  one  another  or 
have  a  specified  separation.  One-dimensional  gap  elements  can  be  used  to  connect  any  two 
elements  in  a  mesh.  The  gap  can  then  be  defined  to  have  a  specific  heat  conductance,  electrical 
conductance,  flow  for  pore  pressure  elements,  heat  generation,  or  radiation.  These  elements 
could  potentially  be  used  to  remove  heat  at  the  location  of  a  cooling  coil  in  the  mesh,  and 
distribute  it  to  an  arbitrary  added  element  that  is  not  part  of  the  actual  model.  This  was  completed 
as  shown  in  Figure  3-2. 

Figure  3-2:  Use  of  a  Gap  Element  (Scale:  maximum  134° F,  minimum  95.1°F,  interval  3°F) 


1 

r 

“■  • 

1 

The  model  shown  had  a  DINTER  type  one-dimensional  heat  transfer  element  connecting  the 
location  of  the  cooling  coil  in  the  center  of  the  mesh  to  an  additional  element  separate  from  the 
mesh.  The  gap  had  a  defined  gap  conductance,  gap  clearance,  average  temperature,  and  average 
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mass  flow  rate  per  unit  area.  The  elemental  interface  also  was  defined  to  have  an  element  cross- 
sectional  area  and  an  x,  y,  and  z-direction  cosine.  The  extra  element  was  given  concrete 
properties  and  an  initial  temperature  equal  to  the  input  temperature  of  the  cooling  coil  water. 

It  can  be  seen  in  Figure  3-2  that  the  use  of  the  gap  element  does  remove  heat  from  the  system  and 
transfers  it  to  the  additional  element  connected  by  the  gap.  However,  it  was  determined  that  it 
was  not  possible  to  modify  this  method  to  model  the  actual  cooling  coil  problem.  The  rate  that 
heat  will  be  transferred  out  of  the  concrete  is  again  dependent  on  the  temperature  differential 
between  the  water  and  concrete.  Using  this  gap  element,  the  temperatures  could  not  be  controlled 
because  there  was  no  subroutine  available  except  GAPCON  to  modify  the  gap  conductance. 

Also,  as  the  additional  element  became  warmer,  the  rate  at  which  heat  would  leave  the  concrete 
would  slow.  This  is  true  in  the  real  life  situation  also,  but  could  not  be  accurately  modeled  in  this 
case.  At  some  point  it  could  even  be  assumed  that  the  additional  element  could  heat  up  to  the 
temperature  of  the  concrete  and  no  longer  work  as  a  cooling  coil.  In  addition,  the  ease  of  use  is 
not  optimal  because  this  option  requires  the  addition  of  several  additional  elements  to  the  system, 
then  the  connection  of  all  the  potential  cooling  coil  nodes  to  the  additional  elements.  It  would 
need  to  be  determined  whether  many  nodes  could  be  gap  connected  to  the  same  element  or  if  the 
rate  of  heat  transfer  would  be  affected  by  the  number  of  gaps  attached.  In  light  of  these  potential 
problems  and  the  apparent  difficulty  of  accurately  modeling  the  actual  problem  using  this 
method,  a  different  method  was  desired. 

3.4  Boundary  Nodes 

Boundary  conditions  can  be  used  to  specify  the  values  of  all  basic  solution  variables  at  nodes 
including  displacements,  rotations,  pore  pressures,  temperatures,  electrical  potentials,  normalized 
concentrations,  or  acoustic  pressures.  They  can  be  assigned  as  ‘model’  input  data  to  define  zero¬ 
valued  boundary  conditions  or  ‘history’  input  data  to  add,  modify,  or  remove  zero-valued  or  non¬ 
zero  boundary  conditions.  Subroutine  ‘DISP’  is  available  to  define  or  redefine  any  magnitudes  of 
prescribed  boundary  conditions. 

Boundary  conditions  can  be  used  in  heat  transfer  analyses  to  specify  boundary  temperatures 
during  the  analysis.  If  these  boundary  temperatures  are  included  at  nodes  within  the  mesh,  they 
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will  act  as  heat  sources  or  sinks  depending  on  the  specific  conditions  of  the  problem.  The 
BOUNDARY  option  was  included  in  the  preliminary  mesh  to  test  the  affects  of  the  boundary 
temperature  as  a  heat  sink.  The  result  is  shown  in  Figure  3-3. 

Figure  3-3:  Preliminary  Mesh  with  BOUNDARY  Condition 

(Scale:  maximum  89.5° F,  minimum  52.1°F,  interval  2.88° F) 


This  plot  demonstrates  the  effect  of  a  heat  sink  after  720  days.  Again,  the  ambient  temperature  is 
applied  to  the  top  and  bottom  boundaries  and  a  symmetry  boundary  to  the  two  sides.  It  is  clear 
from  this  plot  that  the  node  in  the  center  of  the  mesh  with  the  temperature  boundary  condition  is 
acting  to  cool  the  system.  This  method  was  eventually  chosen  to  model  the  coils  and  will  be 
discussed  at  length  in  Chapters  4,  5,  and  6. 

3.5  Forced  Convection  Through  a  Mesh 

Forced  convection  of  a  fluid  through  the  mesh  can  be  modeled  in  ABAQUS.  Forced 
convection/diffusion  heat  transfer  elements  are  used  and  the  velocity  of  the  fluid  moving  through 
the  mesh  is  prescribed  using  the  MASS  FLOW  RATE  option.  This  type  of  analysis  is  intended 
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for  use  in  thermal  problems  involving  heat  transfer  in  a  flowing  fluid  so  that  heat  is  transported 
(convected)  by  the  velocity  of  the  fluid  and  additionally  is  diffused  by  conduction  through  the 
fluid  and  its  surroundings.  The  mass  flow  rate  of  the  fluid  through  the  mesh  affects  the 
conduction  between  the  fluid  and  the  adjacent  forced  convection/diffusion  heat  transfer  elements. 

3.5.1  Peclet  and  Courant  Numbers 


The  numerical  solution  of  the  transient  heat  transfer  equation  including  convection  becomes  more 
complicated  as  convection  dominates  diffusion.  The  Peclet  number  (f),  shown  in  Eq.  3-1,  gives  a 
representation  of  convection  dominance  over  diffusion.  The  Peclet  number  should  be  kept  below 
1000. 


r=\v\M^ 


(3-1) 


Where  |v|  is  the  magnitude  of  the  velocity  vector  and  Al  is  the  characteristic  length  of  the  element 
in  the  direction  of  flow.  If  y=0,  this  implies  no  convection  because  there  is  no  fluid  velocity.  As 
Y->oo,  the  problem  becomes  purely  convective  because  there  is  no  time  for  diffusion.  Petrov- 
Galerkin  finite  element  are  used  in  ABAQUS  to  model  systems  with  high  Peclet  numbers.  These 
elements  control  numerical  diffusion  and  dispersion  by  using  non-symmetric,  upwinded 
weighting  functions  and  stabilize  results. 


Convection/diffusion  elements  with  numerical  dispersion  control  have  a  numerical  stability  limit 
on  the  allowable  time  increment.  This  requirement  is  defined  by  the  Courant  number.  The 
Courant  number  is: 


C  = 


At 

V  — 

Al 


(3-2) 


Where  v  is  the  velocity  of  the  fluid,  At  is  the  time  increment  specified  in  the  analysis,  and  Al  is  the 
characteristic  element  length.  The  Courant  number  measures  how  quickly  energy  can  be 
convected  across  and  element  compared  to  the  time  increment.  If  convective/diffusive  elements 
are  used  in  ABAQUS,  they  cannot  provide  accurate  transient  solutions  for  C  >  1.  If  C  >1, 
energy  can  convect  across  more  than  a  single  element  in  a  time  increment. 
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The  Portuguese  Dam  problem  has  a  flow  rate  through  the  cooling  coil  of  4  gal/min  and  a  0.887” 
inner  pipe  diameter  as  mentioned  earlier.  This  corresponds  to  a  Peclet  number  of  1392  using  two 
foot,  eight-noded  elements  that  have  a  characteristic  length  of  one  foot  (node-to-node  distance). 
However,  long  before  the  Peclet  number  reaches  the  1000  limit  recommended  by  ABAQUS,  the 
Courant  number  exceeds  one  and  the  analysis  has  become  unstable.  Since  the  Courant  number 
limit  is  exceed  at  a  much  lower  threshold  than  the  Peclet  number  for  the  mesh  used  for  the  dam 
problems,  it  is  the  limiting  criteria  for  stability  and  is  discussed  in  greater  length  below. 

The  Courant  number  for  the  Portuguese  Dam  problem  is  45,000  for  the  smallest  time  increment 
used  (1/4  days).  It  increases  for  all  the  larger  time  increments  used  during  the  analysis.  As  stated 
above,  the  Courant  number  should  be  less  than  one.  Having  a  Courant  number  of  45,000  is  a 
huge  variation  from  one,  and  completely  destroys  the  stability  of  the  analysis.  Typically,  the 
water  flowing  through  the  mesh  should  pick  up  heat  from  the  mesh.  Therefore,  the  temperature 
of  the  concrete  surrounding  the  coil  should  gradually  decrease,  and  the  water  in  the  coil  should 
increase  in  temperature.  When  the  solution  breaks  down,  oscillations  are  seen  in  the  temperature 
solutions,  and  the  temperatures  are  not  accurate.  A  3D  analysis  was  completed  on  the  mesh 
shown  in  Figure  3-4.  One  dimensional  convection  diffusion  elements  were  placed  in  an  ‘L’  shape 
in  the  model  to  represent  the  flowing  fluid  through  the  mesh.  The  flowing  fluid  is  cooling  the 
mesh,  as  shown  by  the  temperature  contours  present  in  the  mesh  in  Figure  3-4. 

Figure  3-4:  3D  Mesh  with  Convection/Diffusion  Elements  to  Model  Flowing  Water 
(Scale:  maximum  75° F,  minimum  5(fF,  interval  1.92°F) 
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However,  a  closer  analysis  of  the  cooling  coil  water  temperatures  in  Figure  3-5  shows  the 
breakdown  in  the  solution  to  oscillating,  inaccurate  temperatures  of  the  cooling  water  in  the 
solution.  The  straight  line  at  50°F  on  the  graph  is  the  input  temperature  of  the  water.  The  two 
oscillating  curves  represent  temperatures  of  the  flowing  fluid  at  the  bend  and  exit  point 
respectively  along  the  flow  length. 

Figure  3-5:  Oscillations  in  Destabilized  Solution  due  to  High  Courant  Number 


To  analyze  the  affects  of  elevated  Courant  numbers  on  the  accuracy  of  the  analysis,  the  point  of 
breakdown  in  the  analysis  was  located.  To  do  this,  parameters  of  the  Portuguese  Dam  problem 
were  kept  constant  except  for  the  flow  rate  of  the  water  using  the  preliminary  3D  mesh  and  a 
straight  16’  cooling  coil  as  shown  in  Figure  3-6.  Figure  3-6  shows  one-half  the  model;  the  top  lift 
is  removed  in  order  to  view  the  cooling  coil  placed  between  the  two  lifts.  The  initial  temperature 
of  the  concrete  is  75°F,  the  inlet  temperature  of  the  water  is  SOT,  and  no  heat  generation  is 
included  in  the  concrete. 
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Figure  3-6: 


Mesh  to  Analyze  Courant  Number  Variation  With  16’  Cooling  Coil 
(Scale:  maximum  74.7° F,  minimum  50°F,  interval  2.9F) 


The  mass  flow  rate  corresponding  to  4  gallons/minute  and  a  0.887”  diameter  pipe  is  248,000 
slugs/day/ft^.  This  corresponds  to  the  Courant  number  of  45,000  for  a  Va  day  increment 
mentioned  earlier  for  the  Portuguese  Dam  problem.  When  the  Courant  number  reaches  one,  the 
solution  does  not  immediately  break  down  and  begin  oscillating.  The  solution  begins  to  show 
signs  of  stability  problems  with  Courant  numbers  higher  than  1.6,  and  get  progressively  worse  as 
shown  in  the  figure  series  to  follow.  Figure  3-7  and  3-8  show  stable  solutions  with  Courant 
numbers  of  0.9  and  1.6,  and  mass  flow  rates  of  13.968  slugs/day/ft^  and  24.832  slugs/day/ft^ 
respectively  for  a  five  day  time  period.  The  temperature  plotted  is  the  temperature  of  the  cooling 
water  at  the  cooling  pipe  exit  point  as  shown  in  Figure  3-6.  The  inlet  water  temperature  is  50°F. 


Figure  3-7:  Courant  Number  of  0.9;  Mass  Flow  Rate  of  13.968  slugs/day/ff 


TOTAL  TIME 


Figure  3-8:  Courant  Number  of  1.6;  Mass  Flow  Rate  of 24.832  slugs/day/fi^ 
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The  two  solutions  shown  above  are  stable,  even  though  the  Courant  number  is  1.6  in  Figure  3-8. 
This  can  be  interpreted  to  mean  that  a  cutoff  Courant  number  is  a  safe  guideline  to  follow, 
however  stability  in  the  solution  may  occur  for  Courant  numbers  slightly  higher  than  1.  As  the 
Courant  number  continues  to  increase  with  higher  mass  flow  rates,  instability  is  first  seen  in  the 
descending  portion  of  the  temperature  curve  as  shown  in  Figure  3-9  for  a  Courant  number  of  4.0 
and  mass  flow  rate  of  62.08  slugs/day/ft^. 

Figure  3-9:  Courant  Number  of  4.0;  Mass  Flow  Rate  of 62.06  slugs/day/ft^ 
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The  instability  is  viewed  as  a  loss  of  the  smoothness  on  the  first  portion  of  the  curve  prior  to  the 
peak  and  in  the  descending  portion  of  the  curve  as  slight  oscillations  begin  to  form.  More  serious 
degradation  of  the  solution  is  viewed  by  Courant  number  19.33  corresponding  to  a  mass  flow  rate 
of  300  slugs/day/ft^  as  shown  in  Figure  3-10. 


Figure  3-10:  Courant  Number  of  19.33;  Mass  Flow  Rate  of 300  slugs/day/ft^ 
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TOTAL  TIME 

Almost  total  collapse  of  the  solution  is  viewed  as  location  of  the  peak  temperature  is  severely 
compromised  and  oscillations  are  becoming  more  apparent.  As  the  Courant  number  continues  to 
increase,  complete  collapse  of  the  solution  with  full  oscillation  of  the  temperature  occurs  as 
demonstrated  previously  in  Figure  3-5. 

3.5.2  Limitations  of  Forced  Convection  Through  Mesh 

Several  things  are  of  note  at  this  point.  The  first  is  that  for  the  analysis  shown  in  Figure  3-10,  the 
Peclet  number  is  only  1.2.  This  is  well  below  the  1000  limit  recommended  by  ABAQUS  for 
solution  accuracy.  This  demonstrates  that  in  these  cases,  stability  due  to  the  Courant  number  is 
by  far  the  limiting  case  in  the  analysis.  It  is  also  important  to  recognize  that  by  a  Courant  number 
of  19.33  as  shown  in  Figure  3-7,  the  solution  has  degraded  to  the  point  that  all  accuracy  is  lost. 
This  corresponds  to  a  mass  flow  rate  of  300  slugs/day/ftl  As  stated  earlier,  the  flow  rate  of  the 
Portuguese  Dam  problem  corresponds  to  a  mass  flow  rate  of  248,000  slugs/day/ft^.  Based  on  this 
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comparison,  it  is  not  possible  that  the  Portuguese  Dam  problem  can  be  accurately  modeled  using 
this  method.  One  alternative  considered  was  to  attempt  to  use  similitude  within  ABAQUS  to 
solve  the  Portuguese  Dam  problem  by  decreasing  the  flow  rate  and  modifying  other  parameters, 
such  as  the  conductivity,  density,  or  specific  heat  of  the  concrete  or  water,  in  order  to  obtain  the 
correct  solution.  Several  problems  were  encountered  with  this  method.  The  first  was  that  when 
the  parameters  of  the  water  and  concrete  were  modified,  the  shape  of  the  temperature  curve 
varied  from  the  shape  demonstrated  in  the  above  figures.  This  leads  to  the  conclusion  that  the 
solution  can  no  longer  be  trusted  for  accuracy.  This  type  of  problem  is  demonstrated  in  Figure  3- 
11.  The  problem  modeled  was  the  same  as  the  problems  modeled  above,  however  the  specific 
heat  and  density  of  water  were  changed  to  accommodate  the  decreased  flow  rate  corresponding  to 
a  Courant  number  of  0.4  for  stability.  It  was  a  possibility  that  the  solution  to  this  problem  would 
also  be  the  solution  to  the  actual  problem. 

Figure  3-11;  Courant  Number  of  0.4;  Density  and  Specific  Heat  of  Water  Modified 
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It  is  clear  from  this  plot  heat  is  not  being  absorbed  by  the  cooling  water  in  the  same  manner  as  the 
earlier  plots  leading  to  questions  regarding  the  solution’s  validity.  A  second  problem  encountered 
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was  that  even  if  the  solution  could  somehow  be  obtained  through  modification  of  different 
parameters,  it  would  only  hold  for  a  particular  increment  size.  This  was  not  considered  adequate 
because  the  time  increment  needs  to  be  changed  throughout  the  analysis  or  the  computational 
effort  would  be  too  great. 

3.6  Summary 

Upon  analying  several  possible  methods  for  modeling  the  embedded  cooling  coils,  it  was 
determined  that  the  use  of  boudnary  temperature  nodes  would  be  most  effective.  This  will  be 
described  in  great  detail  in  the  following  chapters.  The  forced  convection  through  a  mesh 
capabilities  of  ABAQUS  also  demonstrated  promise.  Although  this  method  cannot  be  accurately 
utilized  for  models  of  the  magnitude  and  flow  rates  presented  here,  there  is  the  possibility  for 
future  work  in  this  area  to  verify  the  Y-curves. 
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Chapter  4 

Boundary  Temperature  Model 


The  use  of  a  boundary  condition  to  control  the  temperature  at  locations  of  cooling  coils  to  act  like 
a  heat  sink  is  discussed  briefly  in  section  3.4.  However,  to  accurately  use  boundary  temperatures 
as  cooling  coils  heat  sinks,  it  is  imperative  that  the  cooling  coil  boundary  temperature  be  accurate. 
As  already  discussed,  as  the  water  flows  through  the  mesh  it  absorbs  heat  from  the  mesh.  The 
amount  of  heat  removed  from  the  mesh  by  the  water,  and  the  corresponding  increase  in  the  water 
temperature,  is  partially  dependent  on  the  differential  temperature  between  the  concrete  and  the 
water.  It  is  simply  not  enough  to  keep  the  boundary  temperatures  at  the  inlet  temperature  because 
this  would  be  assuming  greater  cooling  than  what  actually  occurs  because  the  cooling  water 
temperature  will  increase  as  the  concrete  temperature  increases.  To  determine  the  water 
temperature,  the  theory  developed  during  the  Boulder  Canyon  Dam  studies  presented  in  Chapter 
1  is  used. 

4.1  Description  of  Boundary  Temperature  Model 
Process 

The  Y-curves  created  during  the  Boulder  Canyon  Dam  analysis  were  used  to  determine  the  water 
temperature.  As  discussed  above,  the  Y-curves  are  used  to  find  the  exit  temperature  of  the  water 
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when  the  two  parameters  described  in  Eq.  2-6  and  2-7  are  known.  These  parameters  are  listed 
again  below: 


KL 

h}t 


where: 

K  =  conductivity  of  concrete 
L  =  length  measure  along  the  cooling  pipe 
c,v  =  specific  heat  of  water  (or  other  cooling  fluid) 
=  density  of  water  (or  other  cooling  fluid) 

Q  -  cooling  fluid  flow  rate 
h  =  diffusivity  of  concrete 
t  =  time  since  cooling  commenced 
D  =  diameter  of  the  cooled  cylinder 


(4-1) 


(4-2) 


As  the  water  absorbs  heat  from  the  concrete,  it  increases  in  temperature.  Therefore,  the  concrete 
at  the  inlet  end  of  the  pipe  is  cooled  more  than  the  concrete  near  the  exit  because  the  water  is 
cooler  at  the  inlet  than  the  outlet.  However,  the  Bureau  of  Reclamation  accounts  for  this  by 
reversing  the  flow  of  water  every  twelve  hours.  Therefore  it  can  be  assumed  that  the  rate  of 
cooling  of  the  concrete  at  all  locations  along  the  length  of  the  pipe  is  relatively  equal.  The  Y- 
curves  can  be  used  to  determine  the  exit  temperature  of  water  since  cooling  has  commenced  and 
this  value  can  be  averaged  with  the  inlet  water  temperature  to  get  an  average  cooling  water 
temperature  for  the  pipe. 


4.1.1  Incremental  Analysis  Procedure 


An  ABAQUS  analysis  is  divided  into  steps  that  are  divided  into  smaller  time  increments.  At  the 
beginning  of  each  increment,  it  is  possible  to  change  a  variety  of  conditions  including  the 
boundary  temperatures.  If  the  boundary  temperatures  are  controlled  by  a  subroutine,  this 
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subroutine  is  called  at  the  beginning  of  each  increment  to  recalculate  the  boundary  temperature 
specified  for  the  corresponding  increment.  As  can  be  seen  from  Eq.  4-1  and  4-2,  if  the  Y-curves 
are  going  to  be  used  to  calculate  the  exit  temperature  of  the  water,  a  variety  of  parameters  need  to 
be  specified.  The  flow  rate  (0),  concrete  conductivity  (K),  specific  heat  of  water  (Ch,),  density  of 
water  (p),  and  concrete  diffiisivity  (h)  are  constant  for  a  particular  problem.  The  length  of  the 
cooled  cylinder  (L)  is  measured  as  the  length  of  coil  in  a  particular  lift.  This  value  and  the 
diameter  of  the  cooled  cylinder  (D)  will  change  as  the  spacing  of  the  coils  and/or  the  lift  size 
changes.  The  time  (t)  is  the  total  time  since  the  start  of  cooling. 

As  stated  earlier,  the  Y-curves  do  not  account  for  heat  generation.  This  becomes  a  problem  when 
trying  to  find  the  exit  temperature.  Once  the  Y-value  is  determined,  Eq.  4-3  is  used  to  determine 
the  exit  temperature. 

ExitWaterTemp.  =  Y  * ( Intial ConcreteTemp.-  IntialWaterTemp.)+ InitialWaterTemp. 

(4-3) 

The  initial  (input)  temperature  of  the  water  is  known  and  the  initial  temperature  of  concrete  is 
typically  the  placement  temperature.  Then,  using  the  time  since  the  start  of  the  analysis,  the  exit 
temperature  of  the  water  at  any  time  after  the  start  of  the  analysis  can  be  found.  However,  when 
heat  generation  is  occurring  in  the  concrete,  this  does  not  give  an  accurate  result  because  heat 
generation  was  not  included  in  the  derivation  of  the  Y-curves. 

One  way  to  deal  with  heat  generation  and  the  Boulder  Canyon  theory  is  to  solve  the  problem  in  a 
series  of  small  increments,  where  each  increment  is  its  own  problem.  The  values  at  the  end  of 
one  increment  become  the  initial  conditions  for  the  next.  Therefore,  at  the  end  of  each  increment 
the  average  concrete  temperature  is  determined.  This  average  concrete  temperature  is  used  as  the 
placement  (initial)  temperature  of  concrete  for  the  next  increment.  This  solution  should  prove  to 
be  relatively  accurate,  except  that  the  X,  Y,  and  Z  Boulder  Canyon  Theory  developed  in  Chapter 
1  was  developed  to  use  the  time  since  the  beginning  of  cooling  at  the  time  (r).  The  theory  was  not 
designed  for  an  incremental  solution,  but  can  still  be  used  reliably.  The  use  of  an  incremental 
solution  is  necessary  to  account  for  the  effects  of  heat  generation.  A  flow  chart  depicting  the 
method  followed  for  the  solution  is  shown  in  Figure  4-1. 


Figure  4-1:  Incremental  Procedure 
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Increment  1 

Boundary  Temperature  =  (input  water  temperature  +  exit  water  temperature)/2 
Concrete  temperature  at  end  of  increment 

Increment  2 

Used  as  concrete  placement 
temperature.  New  exit  water 
temperature  calculated.  Average  water 
temperature  applied  as  new  boundary 
temperature  over  this  increment. 


Use  of  the  boundary  temperature  at  cooling  coil  locations  is  implemented  efficiently.  It  requires 
the  use  of  two  subroutines.  One  is  DISP  and  the  other  is  URDFIL.  DISP  is  called  at  the 
beginning  of  each  increment.  The  user  inputs  certain  constants,  such  as  the  placement 
temperature  of  the  concrete,  the  input  temperature  of  the  water,  and  the  Y-values  needed.  DISP 
calculates  the  exit  temperature  of  water  for  that  particular  increment  based  on  the  appropriate  Y- 
value  and  the  current  temperature  of  concrete.  After  finding  the  exit  temperature  of  the  water, 
DISP  averages  the  exit  and  inlet  temperatures  and  produces  the  average  temperature  to  be  applied 
as  the  boundary  temperature  for  that  increment.  Then  this  process  is  repeated  at  the  beginning  of 
the  next  increment.  The  subroutine  URDFIL  allows  the  user  access  to  the  ABAQUS  results  file 
during  the  analysis.  The  concrete  node  temperatures  are  stored  in  the  results  file,  along  with  all 
other  analysis  data.  The  results  file  must  be  accessed  by  URDFIL  in  order  to  obtain  the  average 
concrete  node  temperatures  at  the  end  of  the  increment  and  the  length  of  the  time  increment. 

The  necessary  Y-values  are  based  on  the  different  lengths  of  coil  and  the  time  increments  used  in 
the  analysis.  The  user  must  enter  all  the  Y-values  necessary  throughout  the  analysis.  The  Y- 
values  will  change  with  the  changing  increment  times  and  also  as  the  size  of  the  lifts,  and 
therefore  the  coils,  change.  More  detail  on  finding  the  appropriate  Y-values  is  given  in  Chapter 
6.  The  current  temperature  of  concrete  is  taken  as  an  average  of  two  comer  nodes  at 
approximately  33.3%  of  the  diameter  of  the  cooled  cylinder  above  and  below  the  cooling  coils. 
All  cooling  coil  nodes  in  a  particular  lift  are  assigned  the  same  average  temperature  concrete 
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nodes  because  the  temperature  in  the  concrete  should  be  relatively  constant  except  for  edge 
effects.  Again,  more  detail  on  choosing  the  appropriate  average  concrete  temperature  nodes  is 
given  in  Chapter  6. 


4.1.2  Average  Concrete  Temperature  Nodes 


As  demonstrated  above,  it  is  important  to  accurately  determine  the  concrete  temperature  at  the 
end  of  one  increment  to  use  as  the  placement  temperature  at  the  start  of  the  next  increment. 
Therefore,  selection  of  the  appropriate  average  concrete  nodes  impacts  the  exit  temperature  of 
water  calculated  from  the  placement  temperature  of  the  concrete  as  shown  in  Eq.  4-1.  Since  the 
placement  temperature  of  concrete  becomes  the  average  concrete  temperature  at  the  end  of  the 
previous  increment,  incorrect  selection  of  this  node  could  lead  to  inaccurate  results.  The  X- 
values  developed  during  Boulder  Canyon  are  used  to  determine  the  mean  concrete  cylinder 
temperature.  An  ABAQUS  3D  analysis  was  completed  for  a  straight  280’  foot  coil  and  the 
temperatures  at  the  nodes  surrounding  the  cooling  coil  were  plotted. 

The  analysis  was  completed  without  including  heat  generation  because  the  X-charts  do  not 
account  for  heat  generation.  As  shown  in  Eq.  4-1,  the  exit  water  temperature  of  water  is 
dependent  on  the  average  concrete  temperature  at  the  end  of  the  increment.  Therefore,  several 
nodes  were  used  to  calculate  the  average  concrete  temperature.  The  nodes  were  comer  nodes  and 
had  increasing  radial  distances  from  the  cooling  coil.  The  model  is  similar  to  the  one  shown  in 
Figure  3-3,  except  that  the  coil  is  280’  long  in  the  ‘1’  direction,  instead  of  16’  long  as  shown  in 
Figure  3-3.  A  drawing  of  the  node  numbers  is  shown  in  Figure  4-2. 


Figure  4-2:  Node  Numbers  for  280’  3D  Mesh 
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The  is  the  portion  of  the  node  number  related  to  the  location  of  the  node  in  the  ‘  1’  direction. 
For  example,  the  nodes  of  the  very  left  end  of  the  280’  model  would  have  a  of  ‘000’  (ie. 
312000).  The  nodes  at  the  halfway  point  are  located  140’  down  the  coil  and  have  a  value  of 
‘140’  (ie.  312140).  The  nodes  numbered  in  Figure  3-10  are  corner  nodes.  Between  each  two 
nodes  shown  is  a  midnode.  This  midnode  has  a  number  designation  between  that  of  the  two 
corner  nodes  (ie.  the  midnode  between  corner  nodes  310140  and  312140  is  31 1 140).  The  280’ 
long  cooling  coil  was  placed  in  the  plane  of  node  320140  which  is  the  plane  of  interface  between 
lift  1  and  lift  2.  The  first  model  results  are  shown  in  Figure  4-3.  This  model  used  nodes  321 140 
and  319140  as  the  two  average  concrete  nodes. 


Figure  4-3:  280’  3D  Model  Using  Nodes  1  ’  Away  As  Concrete  Temperature 
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The  concrete  temperatures  at  the  corresponding  nodes  above  and  below  the  coil  are  equal  to  one 
another.  Therefore,  Figure  4-3  shows  only  the  temperatures  in  Lift  2.  The  lowest  curve,  321 140, 
shows  the  node  temperatures  one  foot  from  the  coil.  Each  line  thereafter  shows  the  concrete 
temperature  at  progressively  increasing  distances  in  one-foot  increments.  Therefore,  325140  is 
the  furthest  node  from  the  cooling  coil  plotted,  at  five  feet  away.  The  line  labeled  ‘X-curves’  is 
the  result  obtained  when  the  result  is  calculated  by  hand  using  the  Boulder  Canyon  X-curves  to 
find  the  mean  concrete  temperature.  Figure  4-3  demonstrates  that  a  node  one  foot  away  is  a  poor 
choice  for  the  average  concrete  temperature  because  the  curve  321140  does  not  match  up  with  the 
‘X-curve’  mean  temperature.  The  result  when  nodes  322140  and  318140  are  used  as  the  two 
average  concrete  nodes  is  shown  in  Figure  4-4. 
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Figure  4-4:  280’  3D  Model  Using  Nodes  2’  Away  As  Concrete  Temperature 


Concrete  Temperatures  at  Nodes 

Average  Concrete  Temperature  Obtained  from  Nodes  2’  from  Cooling  Coils 


x-curves 
-*-321140 
-A- 322140 
-5<- 323140 
-5ie- 324140 
-•-325140 


This  graph  demonstrates  that  using  average  concrete  nodes  2’  away  from  the  cooling  coil  is  also 
not  the  best  choice  because  line  322140  does  not  match  the  ‘X-curves’  mean  temperature  line. 
Figures  4-5,  4-6,  and  4-7  show  the  results  using  average  concrete  nodes  3’,  4’,  and  5’  away 
respectively. 


Figure  4-5:  280’  3D  Model  Using  Nodes  3’  Away  As  Concrete  Temperature 


Concrete  Temperatures  at  Nodes 

Average  Concrete  Temperature  Obtained  from  Nodes  3’  from  Cooling  Coil 
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Figure  4-6:  280’  3D  Model  Using  Nodes  4’  Away  As  Concrete  Temperature 


Figure  4-7:  280’  3D  Model  Using  Nodes  5’  Away  Concrete  Temperature 


From  these  three  graphs,  it  can  be  determined  that  the  mean  concrete  temperature  given  by  the 
‘X-curve’  line  is  best  represented  by  the  concrete  temperatures  3’  and  4’  from  the  cooling  coil. 
This  is  demonstrated  by  the  curves  because  the  3’  and  4’  node  spacings  in  their  respective  graphs 
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most  closely  correspond  to  the  X-curves  temperature  value.  The  spacing  of  the  coil  in  this  model 
is  lO’VxlO’H.  Using  Table  One,  the  diameter  of  the  cooled  cylinder  corresponding  to  this 
spacing  is  11.28’.  Thus,  3’  and  4’  radial  distances  correspond  to  26.6%  and  35.5%  of  the 
diameter  respectively.  This  leads  to  the  recommendation  to  use  the  closest  comer  node  to  33.3% 
of  the  cooled  cylinder  diameter.  It  is  also  important  to  note  that  variations  in  the  choice  of  the 
node  only  marginally  impact  the  resulting  average  boundary  temperature.  This  can  be 
demonstrated  by  the  following  two  graphs.  The  difference  between  using  a  concrete  node  4’ 
away  as  the  average,  which  is  quite  accurate,  and  a  concrete  node  1’  away  as  the  average,  which 
is  not  accurate,  is  inconsequential  as  shown  in  Figures  4-8  and  4-9. 


Figure  4-8:  Average  Boundary  Temperature  Using  Nodes  4  ’  Away  As  Concrete  Temperature 
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This  graph  displays  several  lines  that  are  not  smooth,  but  instead  have  steps.  These  steps  will  be 
discussed  in  more  detail  below.  The  first  line  plotted  is  the  inlet  water  temperature  (which 
remains  constant).  The  exit  water  temperature  is  determined  by  hand  using  the  X  and  Y-curves. 
The  exit  water  temperature  is  averaged  with  the  input  water  temperature  to  obtain  the  average 
water  temperature-curves.  The  ‘avg  water-abaqus’  line  displays  the  ABAQUS  result  obtained 
using  the  average  concrete  temperature  from  the  chosen  nodes.  In  the  case  of  Figure  4-8,  they  are 
nodes  4’  away  from  the  coil.  As  described  above,  the  4’  radial  nodes  are  good  approximations 
for  the  average  temperature  of  concrete.  This  is  observed  by  the  proximity  of  the  ‘Avg  Water- 
curves’  and  ‘Avg  Water-Abaqus’  lines.  This  effect  is  not  seen  when  nodes  1’  away  are  used  for 
the  average  concrete  temperature,  as  shown  in  Figure  4-9. 
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Figure  4-9:  Average  Boundary  Temperature  Using  Nodes  1  ’  Away  As  Concrete  Temperature 


However,  even  though  the  ‘Avg  Water  -curves’  and  Avg  Water  -  Abaqus’  lines  for  the  1’ 
average  concrete  nodes  solution  do  not  line  up  nearly  as  well  as  in  Figure  4-8,  the  difference  is 
about  1°F.  Figure  4-3  shows  that  the  1’  concrete  node  temperatures  are  about  9'’F  from  the 
average  concrete  temperatures  as  determined  by  the  X-curves.  However,  this  only  translates  to 
about  a  TF  difference  in  average  water  temperature.  This  water  temperature  then  impacts  the 
calculated  concrete  temperatures  for  the  next  increment  only  marginally.  This  is  important 
because  the  primary  concern  is  obtaining  accurate  concrete  temperatures.  Therefore,  using  the 
closest  comer  node  to  (33.3%)  is  accurate  enough  for  most  purposes. 

4.1.3  Incremental  vs.  Total  Time  Analysis 

The  last  important  issue  to  be  addressed  with  the  boundary  solution  is  the  steps  seen  in  the  water 
temperature  as  demonstrated  by  Figures  4-8  and  4-9.  These  steps  illustrate  a  very  important  point 
-  that  the  curves  developed  during  the  Boulder  Canyon  Dam  investigation  are  not  meant  to  be 
used  to  solve  an  incremental  analysis.  This  means  that  the  t  value  in  Equation  2-7  is  supposed  to 
be  the  total  time  since  the  start  of  cooling,  as  stated  earlier.  However,  since  the  curves  do  not 
include  heat  generation,  it  is  impossible  to  incorporate  heat  generation  into  the  solution  without 
using  the  incremental  solution  described  above,  where  each  time  period  is  solved  as  its  own  small 
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problem,  and  the  ending  parameters  from  one  increment  are  used  as  the  initial  values  for  the  next. 
When  using  ABAQUS,  these  time  increments  can  be  changed  throughout  the  analysis.  When  the 
time  increment  changes,  the  Y-value  used  changes  because  the  t  value  in  Eq.  2-7  changes.  This  is 
why  the  user  is  required  to  enter  a  new  Y-value  for  the  different  time  increments,  as  mentioned 
earlier.  At  each  Y-value  change,  a  jump  is  seen  in  the  exit  temperature  values  found  using  Eq.  4- 
1 .  These  are  the  jumps  seen  in  Figures  4-8  and  4-9. 

To  analyze  the  error  in  the  solution  when  an  incremental  solution  is  used  instead  of  a  total  time 
solution,  the  problem  without  heat  generation  was  solved  several  different  ways.  Using  a  280’ 
straight  coil,  3D  model  with  10’  x  10’  lift  spacing  and  no  heat  generation  the  problem  was  solved 
using  ABAQUS.  Two  types  of  ABAQUS  analysis  were  completed.  The  first  was  a  total  time 
analysis.  The  Y-values  used  changed  for  every  increment  as  the  total  time  increased.  The  only 
difference  between  solving  the  analysis  using  ABAQUS  in  this  manner  vs.  completing  the 
analysis  by  hand,  is  that  ABAQUS  calculates  the  concrete  temperatures.  Otherwise  the  concrete 
temperatures  would  need  to  be  found  using  the  X  and  Z-charts.  For  the  incremental  analysis,  the 
procedure  was  the  same  as  that  described  above.  The  analysis  was  completed  in  Va  day 
increments  which  is  considered  to  be  the  minimum  increment  specified  by  the  Army  Corp  of 
Engineers  for  NISA  models.  The  average  concrete  temperature  was  based  either  on  nodes  3  feet 
away  (Figure  4-10)  or  nodes  4  feet  away  (Figure  4-11). 


Temperature  (dej 


54 


Figure  4-10:  Total  Time  vs.  Incremental  Analysis:  ABAQUS  Incremental  Temperatures  Based 
on  Average  Concrete  Nodes  Three  Feet  Away 
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Both  Figure  4-10  and  4-1 1  are  essentially  the  same,  again  demonstrating  the  insignificant 
difference  made  by  the  choice  of  the  average  concrete  nodes  as  long  as  they  are  close  to  the 
correct  values  as  predicted  by  the  X-curves.  The  nodes  numbers  321 140,  322140,  323 140, 
324140,  and  325140  signify  concrete  nodes  at  140’  along  the  length  of  the  coil  and  at  1’,  2’,  3’, 

4’,  and  5’  radial  distances  away  from  the  cooling  coil  respectively.  The  two  curves  ‘Avg.  Water- 
ycurves-total’  and  ‘Avg.  Water-Abaqus-total’  are  the  average  of  the  inlet  and  exit  water 
temperatures  from  a  hand  calculation  analysis  using  the  Y-curves  and  the  ABAQUS  result 
respectively.  Only  one  curve  is  visible  because  the  curves  are  the  same  and  plot  on  top  of  one 
another,  which  is  expected.  The  two  curves  ‘Avg.  Water  -  ycharts  -  inc’  and  ‘Avg.  Water  - 
Abaqus  -  inc’  are  the  average  water  temperatures  determined  with  an  incremental  analysis  from  a 
hand  calculation  analysis  using  the  Y-curves  and  the  ABAQUS  results  respectively.  The  slight 
difference  in  these  two  curves  is  due  to  the  difference  in  the  average  concrete  node  chosen  in 
ABAQUS  (‘323140-inc.’  or  ‘324140  -  inc.’)  vs.  the  X-curve  average  concrete  value  ‘Concrete 
Temp  -  xcurves  -  total’.  Again,  it  can  be  seen  the  difference  in  the  average  water  curves  is  slight. 
Of  most  importance  is  the  error  between  the  average  water  temperatures  calculated  using  the 
incremental  or  the  total  time  analysis.  These  values  have  a  maximum  difference  of  1.3°  (Figure 
4-10).  This  is  not  necessarily  insignificant,  however  it  leads  to  very  small  differences  in  the 
concrete  temperatures  which  is  what  is  ultimately  the  most  important.  The  concrete  temperatures 
plotted  at  each  node  for  the  incremental  and  total  time  analysis  plot  right  next  to  one  another,  with 
the  most  discrepancy  seen  at  the  node  closest  to  the  cooling  coil,  321 140.  This  difference  has  a 
maximum  value  of  0.52°.  Therefore,  it  is  concluded  that  the  error  incurred  in  the  concrete 
temperature  results  by  using  incremental  analysis  techniques  instead  of  total  time  analysis  is 
minimal. 

4.2  Summary 

The  boundary  node  method  for  modeling  the  cooling  coils  has  several  advantages,  and  also 
disadvantages.  A  major  advantage  of  using  the  boundary  method  of  modeling  the  cooling  coils  is 
that  it  is  easily  implemented  into  an  existing  model.  No  change  of  the  mesh  is  required,  and  no 
special  mesh,  element  size,  or  node  placement  is  required  to  use  the  procedure.  The  actual 
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mechanics  of  implementing  the  procedure  into  a  model  are  outlined  in  Chapter  6.  However,  the 
implementation  is  not  rigorous,  and  requires  no  subroutine  modification  by  the  user  beyond 
entering  certain  parameters  such  as  Y-values,  placement  concrete  temperature,  and  input  water 
temperature.  Also,  implementation  is  essentially  the  same  in  the  two-dimensional  and  three- 
dimensional  analysis,  reducing  confusion  for  the  user. 

The  two-dimensional  model  would  be  most  likely  used  to  try  a  variety  of  cooling  coil  spacings 
and  designs.  It  is  easy  to  modify  the  cooling  coil  spacings  as  long  as  the  user  has  nodes  placed 
where  they  would  like  to  place  coils.  In  addition,  the  three  dimensional  model  can  accommodate 
any  winding  pattern  for  the  cooling  coils  layout.  Another  important  point  to  mention,  which  will 
be  discussed  again  in  Chapter  6,  is  that  the  model  does  not  require  that  the  cooling  coils  be 
functional  for  the  duration  of  the  analysis.  The  coils  can  be  turned  on  for  any  steps  in  the 
analysis,  then  turned  off  when  the  user  would  like  to  stop  cooling.  They  can  also  be  turned  back 
on  later  in  the  analysis  if  desired  by  the  user. 

The  most  important  disadvantage  to  the  boundary  node  method  of  modeling  the  cooling  coils  is 
the  reliance  on  the  Y-curves.  The  Y-curves  were  developed  based  on  several  assumptions, 
including  an  insulated  cylinder,  a  straight  coil,  and  non-incremental  analysis.  However,  the 
curves  have  been  used  extensively  by  the  Bureau  of  Reclamation  since  the  construction  of  the 
Hoover  Dam,  with  good  results.  In  addition,  as  has  been  established,  the  error  due  to  an 
incremental  analysis  is  minor.  However,  a  method  of  analysis  which  does  not  require  the  use  of 
the  Y-curves  would  be  desirable. 

The  other  disadvantage  of  using  boundary  nodes  to  model  the  coils  is  that  the  coils  must  be 
placed  at  node  locations.  Due  to  this  constraint,  the  user  should  consider  possible  spacing  options 
before  creating  their  mesh  in  order  to  choose  a  node  spacing  which  allows  for  analysis  of  several 
coil  spacing  options. 
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Chapter  5 

Results  of  Boundary  Temperature  Model 


This  chapter  will  demonstrate  the  results  of  the  cooling  coil  model  implemented  in  actual  dam 
systems.  The  chapter  begins  with  simplistic  representative  problems  that  demonstrate  the 
effectiveness  of  the  model,  both  in  2D  and  3D.  The  last  model  presented  is  the  2D  Portuguese 
Dam  model  from  the  Army  Corp  of  Engineers. 

5.1  Preliminary  Results 

First,  a  two-dimensional  model  will  be  presented.  This  model  includes  the  sequential  placement 
of  lifts,  films,  heat  generation,  and  other  ‘real  life’  parameters.  Then,  an  exampleJD  model  will 
be  presented  in  less  detail.  Finally,  two  models  will  be  presented  which  demonstrate  the 
compatibility  of  the  2-dimensional  and  3-dimensional  results. 

5.1.1  Two-Dimensional  Model  Results 


To  demonstrate  the  model  results,  a  2D  model  is  presented  below.  This  model  has  three  lO’V  x 
16’H  lifts  placed  sequentially  in  seven  day  increments.  Film  coefficients  surround  the  boundaries 
of  the  lift  exposed  to  the  air.  The  film  coefficients  are  modeled  with  the  amplitude  ambient 
temperature  for  the  Portuguese  Dam.  The  base  of  the  first  lift  has  a  symmetry  boundary  to 
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represent  the  interface  with  the  foundation,  which  is  very  conservative.  Heat  generation  is 
included  in  all  three  lifts,  with  the  amount  of  heat  generated  dependent  on  the  time  since  the  lift 
has  been  placed.  Heat  generation  is  controlled  by  the  subroutine  HETVAL.  Figure  5-1 
represents  the  element  configuration  of  the  model,  and  indicates  certain  node  numbers  on  a 
vertical  line  through  the  center  of  the  model.  The  temperatures  at  these  node  locations  will  be 
used  to  later  to  plot  temperature  distributions  in  the  model,  so  it  is  important  to  note  their 
placement  in  the  model. 

Figure  5-1:  Two-Dimensional,  three  lift  model 
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This  model  begins  with  the  first  lift  being  placed  on  day  one,  the  second  on  day  7,  and  the  third 
on  day  14.  The  model  was  run  for  421  days.  The  increment  time  varied,  with  smaller  increments 
being  used  for  the  beginning  portion  of  the  analysis  when  heat  generation  was  occuring,  and 
larger  increments  being  used  as  the  temperature  became  more  constant.  The  cooling  coils  were 
controlled  by  the  boundary  command  and  the  subroutines  DISP  and  URDFIL  as  described  in 
Chapter  4.  Three  cooling  coils  were  placed  at  four  feet  intervals  at  the  interfaces  between  the 
lifts.  The  cooling  coils  were  not  activated  until  the  upper  lift  bordering  the  interface  was  placed. 
Figure  5-2  shows  the  temperature  contour  after  1.25  days  when  only  one  lift  is  in  place  and  no 
cooling  coils  have  been  activated.  The  placement  temperature  of  the  concrete  at  day  0  was  75°. 
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Figure  5-2:  Three  Lift,  2D  Analysis  after  1.25  days 

(Scale:  maximum  100°F,  minimum  90.^F,  interval  O.TF) 


Figure  5-2  shows  that  the  maximum  temperature  reached  after  1 .25  days  due  to  heat  generation  is 
100  degrees.  The  coolest  temperature  is  about  91“  at  the  three  boundaries  being  cooled  by  the 
ambient  temperature.  Again,  the  bottom  edge  of  the  lift  (which  would  be  the  interface  with  the 
foundation  in  an  actual  dam)  has  a  symmetry  boundary.  Figure  5-3  shows  the  model  at  time  7.25 
days;  0.25  days  after  lift  two  has  been  placed.  The  cooling  coils  between  lift  one  and  two  were 
activated  when  lift  two  was  placed  at  day  7. 

Even  after  0.25  days,  the  cooling  of  the  concrete  by  the  coils  is  visualized.  The  maximum 
temperature  in  lift  one  is  1 17°.  The  temperature  in  lift  2  is  roughly  80°.  Lift  2  had  a  placement 
temperature  of  75°  and  has  been  experiencing  heat  generation  for  only  0.25  days.  The  cooling 
coils  at  the  interface  between  lift  1  and  lift  2  shown  here  are  at  a  temperature  of  51.6°.  This 
temperature  is  calculated  as  discussed  in  Chapter  4  by  subroutines  ‘DISP’  and  ‘URDFIL’  using 
an  input  water  temperature  of  40°.  The  average  concrete  nodes  used  are  2240  and  1840  (see 
Figure  28).  Nodes  2240  and  1840  are  at  2’  from  the  cooling  coil.  The  diameter  of  the  cooled 
cylinder  for  lO’V  x  4’H  cooling  coil  spacing  is  7.14’  (Table  1).  Following  the  guidelines  for 
choosing  the  correct  average  concrete  nodes  as  detailed  in  section  4.1.2,  we  have  a  choice 
between  comer  nodes  2’  or  4’  from  the  coil.  This  leads  to  28%  and  56%  of  the  diameter 
respectively.  Therefore,  nodes  2’  away  were  the  best  choice.  They  are  chosen  in  the  center  of  the 
lift  to  avoid  end  effects.  The  same  cooling  coil  temperature  is  applied  to  all  three  coils  since  in 
the  actual  problem  this  would  be  the  same  coil  twisting  through  the  lift.  The  exterior  boundaries 


led  by  the  ambient  temperature.  The  contours  caused  by  cooling  due  to  the  coils  can 
en  at  this  point.  They  more  pronounced  in  Figure  5-4,  after  13  days. 


e  5-3:  Three  Lift,  2D  Analysis  after  7.25  days 

(Scale:  maximum  1  IT F,  minimum  5LTF,  interval  5°F} 


e  5-4:  Three  lift,  2D  Analysis  after  13.0  days 
maximum  120° F,  minimum  53.4° F,  interval  5.1° F) 
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In  Figure  5-4,  the  effect  of  the  cooling  coils  can  be  clearly  seen.  Their  action  is  to  cool  the  center 
portion  of  the  model.  The  highest  temperature  experienced  is  close  to  the  symmetry  boundary,  at 
120°.  The  maximum  temperature  experienced  in  lift  2  is  about  1 10°.  The  area  close  to  the 
cooling  coils  is  85  to  90°F.  The  cooling  coil  temperature  at  this  point  is  53.4°,  It  is  almost  2° 
higher  than  the  cooling  coil  temperature  at  7.25  days  because  the  temperature  of  the  average 
concrete  nodes  is  increasing  due  to  heat  generation.  At  day  14,  lift  3  is  placed.  Figure  5-5  shows 
the  model  at  time  14.25  days,  0.25  days  after  lift  3  has  been  placed. 

Figure  5-5:  Three  Lift,  2D  Analysis  after  14.25  days 

(Scale:  maximum  I20°F,  minimum  51.6°F,  interval  5.25°F) 


Figure  5-5  shows  the  model  with  lift  3  in  place.  Lift  3  is  initially  much  cooler  than  lifts  1  and  2 
because  it  has  only  been  generating  heat  for  0.25  days.  The  cooling  coils  are  activated  between 
both  lifts  one  and  two  and  lifts  two  and  three.  The  coils  at  the  interface  between  lift  two  and  three 
use  3240  and  2840  (see  Figure  28)  as  the  two  average  concrete  nodes  for  the  same  reasons 
described  above  for  nodes  2240  and  1840.  The  cooling  coils  at  the  interface  between  lift  1  and  2, 
and  lifts  2  and  3  can  be  at  different  temperatures  because  their  average  concrete  nodes  are  at 
different  temperatures.  The  difference  will  become  less  pronounced  as  the  heat  generation  ceases 
and  the  concrete  temperatures  became  more  evenly  distributed  throughout  the  model.  The 
contours  after  21.0  days  are  shown  in  Figure  5-6. 
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Figure  5-6:  Three  Lift,  2D  Analysis  after  21.0  days 

(Scale:  maximum  116°F,  minimum  5 1  .T F,  interval  4.95°F) 


In  Figure  5-6,  the  cooling  from  the  coils  is  clearly  visualized.  The  difference  in  the  progression 
of  the  cooling  at  the  interface  between  lifts  1  and  2  compared  to  lifts  2  and  3  is  also  clear.  The 
highest  temperatures  experienced  at  this  point  are  at  the  symmetry  boundary  in  lift  1,  and  in  the 
center  of  lift  three.  Lift  two  is  being  cooled  by  cooling  coils  on  either  side.  At  day  41,  the  extent 
of  cooling  in  lift  2  is  even  more  pronounced,  emphasizing  the  effect  that  the  cooling  coils  have  on 
the  model  as  shown  in  Figure  5-7. 

The  center  of  lift  two  in  Figure  5-7  is  at  about  80®  at  this  point.  This  contrasts  with  the  maximum 
temperature  of  99.7°  at  the  symmetry  boundary,  which  is  experiencing  little  cooling  from  either 
the  coils  or  the  ambient  temperature.  The  maximum  temperature  in  lift  3  is  about  91°.  After 
421.0  days,  the  model  has  cooled  to  the  point  that  the  ambient  temperature  is  higher  than  the 
concrete  and  the  external  temperature  is  actually  heating  the  sides  of  the  model.  This  is  seen  in 
Figure  5-8. 
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The  temperature  in  the  center  of  the  model  is  about  67°,  demonstrating  the  extent  of  cooling 
possible  by  the  coils.  The  boundary  temperature  is  about  85°  around  all  three  sides  due  to  the 
ambient  temperature.  The  progression  of  cooling  is  more  clearly  seen  when  the  node 
temperatures  are  plotted  over  time.  Figure  5-9  shows  the  cooling  coils  temperatures  at  the 
interface  of  each  lift  by  plotting  cooling  coil  nodes  2040  and  3040  over  time. 

Figure  5-9:  Three  Lift,  2D  Model :  Cooling  Coil  Temperatures 
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The  node  number  of  the  cooling  coil  is  given  in  the  box  on  the  upper  left  of  Figure  5-9.  The  steps 
viewed  in  the  lower  part  of  both  curves  are  due  to  the  changing  Y-values  discussed  in  4.3.1.  For 
the  first  seven  days,  the  node  3040  does  not  exist  in  the  model  because  lift  3  has  not  been  placed. 
Therefore  ABAQUS  plots  the  75°  placement  temperature  for  node  3040  for  the  first  7  days.  At 
day  7,  lift  2  is  placed,  and  the  node  becomes  active.  For  both  coils,  during  the  first  seven  days  of 
their  placement  they  function  as  normal  concrete  nodes.  They  experience  heat  generation, 
obtaining  a  maximum  temperature  of  93.0°.  After  they  have  been  in  place  for  7  days,  the  lift 
above  them  is  placed,  and  they  are  abruptly  turned  to  cooling  coils  controlled  by  subroutines 
DISP  and  URDFIL  as  described  in  Chapter  4.  After  heat  generation  ends  and  the  concrete 
temperatures  throughout  the  model  become  more  constant,  the  cooling  coil  temperatures  begin  to 
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level  out,  and  are  almost  equal  at  day  140.  The  cooling  coil  temperatures  compared  to  the 
concrete  temperatures  for  lift  two  are  viewed  in  Figure  5-10. 

Figure  5-10:  Three  Lift,  2D  Model:  Lift  2  Temperatures 
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The  node  number  given  in  the  upper  left  comer  of  Figure  5-10  correspond  to  the  node  numbers  in 
Figure  5-1.  The  two  cooling  coils  2040  and  3040  are  the  same  as  shown  in  Figure  5-1  as  well. 
The  peak  concrete  temperatures  are  affected  by  the  node  proximity  to  a  cooling  coil.  The 
maximum  temperature  reached  was  1 12.9°  at  node  2640.  The  oscillations  that  can  be  seen  in  the 
later  part  of  the  curve  (after  day  80)  occur  because  the  overall  concrete  temperature  has  reduced 
to  the  point  that  it  now  oscillates  according  to  the  ambient  temperature.  The  same  plot  was 
created  for  a  model  with  a  constant  Y-value  instead  of  the  correct,  changing  Y-values  which  lead 
to  steps  in  the  cooling  coils  temperature.  This  plot  is  shown  in  Figure  5-11. 
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Figure  5-11:  Three  Lift,  2D  Model:  Lift  2  Temperatures  with  Constant  Y-value 


The  maximum  concrete  temperature  reached  at  node  2640  using  a  constant  Y-value  is  1 13.0°. 

This  differs  from  the  112.9°  temperature  found  when  using  the  correct,  changing  Y-values  in 
Figure  5-10.  This  again  demonstrates  that  the  error  incurred  by  the  use  of  an  incremental  analysis 
with  changing  Y-values  is  minimal.  The  temperature  plot  for  lift  1  (with  the  correct,  changing  Y- 
values)  is  shown  in  Figure  5-12. 

The  maximum  temperature  in  lift  1  is  119.9°,  reached  by  node  1040  which  is  on  the  symmetry 
boundary.  From  this  upper  boundary,  we  see  a  decreasing  temperature  experienced  by  each  node 
as  we  progress  toward  the  cooling  coil  (ie.  1240,  1440,  etc.).  Again,  we  see  the  oscillations  in  the 
later  parts  of  the  curve  due  to  the  ambient  temperature.  The  lift  3  concrete  temperatures  are 
shown  in  Figure  5-13. 
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The  maximum  temperature  reached  in  lift  3  is  1 12.9°  by  node  3640.  However,  in  lift  3,  the 
effects  of  the  ambient  temperature  are  clearly  seen.  Node  4040,  which  is  sitting  on  the  ambient 
boundary  (see  Figure  5-1),  never  reaches  a  maximum  temperature  close  to  those  experienced  by 
the  interior  portions  of  the  lift.  This  occurs  because  the  boundaries  increase  in  temperature  until 
they  reach  the  ambient  temperature,  at  which  point  they  begin  to  follow  the  ambient  temperature 
curve  as  they  are  cooled  (or  heated)  by  the  external  temperature.  This  effect  is  also  seen  at  node 
3840  whose  peak  is  diminished  by  ambient  temperature  cooling.  The  peak  temperature  for  node 
3240  is  strongly  affected  by  cooling  due  to  coil  3040. 

5.1.2  Three-Dimensional  Model  Results 

A  three-dimensional  model  is  demonstrated  here  in  less  detail  to  present  the  ease  of  modeling  a 
variety  of  cooling  coil  patterns  in  three  dimensions.  The  model  presented  is  a  three  dimensional 
10’Vxl6’Hxl0’D  model  with  two  lifts,  both  in  place  at  the  start  of  analysis.  The  model  has 
ambient  temperatures  applied  at  the  right  and  left  faces,  and  symmetry  boundaries  on  the  bottom, 
top,  front,  and  back  faces  of  the  two  lifts.  In  the  figures  shown,  the  upper  lift  is  removed  so  that 
the  temperature  distribution  can  be  easily  seen.  The  model  includes  a  cooling  coil  that  enters  on 
the  back  face,  winds  through  the  model,  and  exits  on  the  front  face.  Heat  generation  was 
included,  with  the  initial  water  temperature  at  50°  and  the  concrete  placement  temperature  at  75°. 

Figure  5-14  shows  a  maximum  temperature  of  106°  after  1.75  days.  The  cooling  coil  water  is 
51.6°  and  cooling  due  to  the  coils  is  apparent.  The  ambient  temperature  has  cooled  the  left  and 
right  faces  to  approximately  91°.  Close  to  the  cooling  coils,  concrete  temperatures- are  closer  to 
76°.  The  contour  plot  after  7  days  is  shown  in  Figure  5-15. 
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temperature  in  the  lift  is  1  IT  at  the  bottom  symmetry  boundary  and  the  cooling  coil  water  is 
51.4°.  More  extensive  cooling  is  seen  at  days  107  and  257  shown  in  Figures  5-16  and  5-17. 

Figure  5-16:  3D  Analysis  After  107.0  Days 

(Scale:  maximum  74.1°F,  minimum  50.3° F,  interval  1.8°F) 


Figure  5-17 :  3D  Analysis  After  257. 0  Days 

(Scale:  maximum  73.4°F,  minimum  50. 2° F,  interval  1.78°F) 
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In  Figure  5-16,  cooling  contours  due  to  the  cooling  coil  are  extending  well  into  the  lift  after  this 
time  has  passed.  The  concrete  has  cooled  to  the  point  that  the  ambient  temperature  is  warming 
the  concrete  surface  at  the  right  and  left  boundaries,  which  are  approximately  at  74°.  Near  the 
cooling  coils  the  temperatures  are  close  to  60°. 

After  257  days  in  Figure  5-17,  the  effect  of  the  cooling  coils  reaches  down  to  the  bottom 
symmetry  boundary.  The  highest  concrete  temperatures  are  at  the  left  and  right  faces  with  the 
ambient  temperature  applied,  at  about  73°.  Close  to  the  cooling  coil,  the  concrete  temperature  is 
about  55°. 

5.1.3  Comparison  of  2D  and  3D  Model  Results 


It  was  desirable  to  compare  two-dimensional  and  three-dimensional  results  for  a  similar  model  to 
ensure  the  validity  of  the  analysis  in  ABAQUS.  This  process  is  described  below. 

2D  Model  Set-up 

The  2-dimensional  mesh  is  designed  as  shown  in  Figure  2-5.  It  includes  two  lO’V  x  16’H  lifts. 
The  lifts  are  both  in  place  at  the  start  of  the  analysis;  they  are  not  placed  sequentially.  The 
boundaries  were  all  insulated,  no  ambient  temperatures  were  applied.  This  was  done  so  that  the 
2D  and  3D  modeling  of  the  cooling  coils  alone  could  be  observed  and  compared.  Three  coils 
placed  at  four-foot  intervals  in  the  center  of  the  mesh  were  modeled.  Heat  generation  was 
included  and  the  placement  temperature  and  the  input  water  temperature  were  50°  and  75° 
respectively.  The  two-dimensional  analysis  after  five  days  is  shown  in  Figure  5-18.  After  five 
days,  the  maximum  concrete  temperature  is  1 14°.  The  cooling  coil  temperature  is  55.3°  and  the 
contours  created  by  the  cooling  coils  are  apparent. 
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3D  Model  Set-up 

The  3D  model  is  exactly  the  same  as  the  2D  model  except  that  there  is  a  depth  dimension  of  10’. 
The  model  includes  heat  generation  and  insulated  boundaries,  with  placement  concrete  and 
cooling  input  water  temperatures  of  75°  and  50°  respectively.  After  five  days,  the  temperature 
contour  is  shown  in  Figure  5-19.  A  slice  was  taken  through  the  model  in  order  to  view  the 
temperature  contours  within  the  concrete  lift  section. 

Figure  5-19:  Three  Dimensional  Analysis  After  Five  Days  (for  comparison) 

(Scale:  maximum  114°F,  minimum  SS.f’F,  interval  4.5°F) 


The  maximum  temperature  in  the  concrete  is  1 14°,  exactly  the  same  as  the  2D  model.  The 
cooling  water  temperature  is  55.3°,  also  consistent  with  the  2D  model. 

Results 

The  following  sequence  demonstrates  the  compatibility  of  the  2D  and  3D  models. 


Figure  5-20:  2D  and  3D  Analysis  After  20.0  Days 

(Scale:  maximum  123°F,  minimum  53.8°F,  interval  5.3‘’F) 


Figures  5-20  and  5-21  both  show  identical  temperature  distributions  for  the  two-dimensional  and 
three-dimensional  models,  verifying  the  validity  of  their  compatibility. 


5.2  Two-Dimensional  Portuguese  Dam  Results 


The  cooling  coil  modeling  procedure  was  implemented  in  existing  two-dimensional  models  of  the 
Portuguese  Dam  provided  by  the  Army  Corp  of  Engineers.  The  foundation  and  first  four  lifts 
were  analyzed.  The  model  has  2.5’  horizontal  elements  in  the  lifts  with  varying  vertical  element 
height.  Three  cooling  coils  were  placed  at  ten  foot  increments  at  the  lift-lift  or  lift-foundation 
boundaries.  The  lifts  vary  in  size  with  horizontal  dimensions  of  40’  and  varying  vertical 
dimensions.  The  ambient  temperature  was  applied  to  all  exposed  faces  of  the  lifts.  Heat 
generation  was  included  and  the  input  water  temperature  was  50°.  Figure  5-22  shows  the  analysis 
after  7  days. 

Figure  5-22:  Portuguese  Dam  At  7.0  Days 

(Scale:  maximum  103°F,  minimum  55“ F,  interval  3.TF) 


From  the  start  of  the  analysis  the  cooling  coils  between  lift  1  and  the  foundation  are  activated. 
After  7  days,  the  maximum  temperature  in  lift  1  is  1 10°.  The  cooling  water  is  at  53.7°.  Contours 
due  to  cooling  by  the  coils  and  the  ambient  temperature  are  visible.  At  7  days,  lift  two  is  placed. 
The  temperature  after  14  days  is  shown  in  Figure  5-23. 
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Figure  5-23:  Portuguese  Dam  At  14.0  Days 

(Scale:  maximum  115°F,  minimum  53.3°F,  interval  4.75°F) 


After  14  days,  the  maximum  temperature  is  experienced  in  lift  2  at  1 15°.  Both  sets  of  cooling 
coils  are  working  to  cool  lift  one,  whose  temperature  varies  from  approximately  93°  at  the 
foundation  to  108°  closer  to  lift  2.  The  ambient  temperature  is  cooling  the  exposed  boundaries. 
At  14  days,  lift  3  is  placed.  At  21  days,  the  temperature  distribution  is  as  shown  in  Figure  5-24. 


Figure  5-24:  Portuguese  Dam  Analysis  At  21.0  Days 

(Scale:  maximum  115°F,  minimum  53° F,  interval  4.75° F) 
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At  21  days,  the  temperatures  of  both  lift  two  and  three  are  elevated,  at  about  1 15°.  The  ambient 
temperature  is  cooling  the  boundaries,  and  the  cooling  coils  are  also  making  an  impact.  Lift  4  is 
placed  at  21  days,  and  after  28  days,  the  temperature  contours  are  as  shown  in  Figure  5-25. 

Figure  5-25 :  Portuguese  Dam  Analysis  After  28.0  Days 

(Scale:  maximum  1 1 5° F,  minimum  52.8° F,  interval  4.78° F) 


At  28  days,  the  highest  temperatures  are  experienced  in  lift  4,  at  1 15°.  Lifts  2  and  3  are  similar  in 
temperature;  about  108°.  Lift  1  is  significantly  cooler  due  to  the  cooling  coils,  and  the  effects  of 
the  ambient  temperature  cooling  are  observed  along  the  boundaries.  As  cooling  continues,  the 
temperature  continues  to  diminish,  as  seen  in  Figure  5-26  at  53  days. 

Figure  5-26:  Portuguese  Dam  Analysis  at  53.0  Days 

(Scale:  maximum  lll°F,  minimum  51.TF,  interval  4.56°F) 
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After  53  days,  the  extent  of  cooling  is  more  extensive.  Lift  1  is  at  approximately  77°,  lifts  2  and  3 
are  about  90°,  and  the  maximum  temperature  is  experienced  at  the  top  of  lift  4  at  1 1 1°.  At  278 
days,  the  cooling  has  advanced  substantially,  as  seen  in  Figure  5-27. 

Figure  5-27:  Portuguese  Dam  Analysis  at  278.0  Days 

(Scale:  maximum  79.8°F,  minimum  SO.S^F,  interval  2.25°F) 


Cooling  at  this  stage  is  quite  extensive.  The  lifts  have  equilibrated  in  temperature,  and  the  central 
portion  near  the  coils  is  about  55°.  The  ambient  temperature  at  the  external  boundaries  is  acting 
to  warm  the  concrete  surfaces  at  this  point. 

5.3  Summary 

The  modeling  of  embedded  cooling  coils  using  boundary  temperature  nodes  can  be  completed 
accurately  and  efficiently  with  good  results  as  presented  in  this  chapter.  The  modeling  procedure 
can  be  modified  to  work  for  any  massive  concrete  problem  and  is  relatively  easy  to  implement, 
which  will  be  discussed  in  Chapter  6.  The  cooling  coil  nodes  effectively  cool  the  massive 
concrete  system  and  errors  due  to  the  use  of  incremental  analysis  and  average  concrete  node 
choice  are  insignificant  as  presented  in  Chapter  5. 
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Chapter  6 

Documentation  for  Cooling  Coil  Modeling 
Using  ABAQUS  Version  5.7  or  5.8 


This  chapter  will  cover  the  use  and  implementation  of  the  commands  and  subroutines  required  to 
model  cooling  coils  in  massive  concrete  systems  using  ABAQUS.  This  chapter  will  discuss  the 
additions  that  need  to  be  made  to  an  ABAQUS  input  file,  the  additional  subroutines  required,  and 
how  the  subroutines  perform  the  necessary  tasks. 

6.1  Basic  Procedure 

The  procedure  developed  allows  the  user  to  implement  a  model  in  ABAQUS  that  analyzes  how  a 
cooling  coil  draws  heat  out  of  a  massive  concrete  system.  It  allows  a  user  to  change  the  spacing 
and  number  of  the  coils  in  the  system  or  a  particular  lift  and  to  turn  the  coils  on  or  off  during 
different  parts  of  the  analysis.  It  also  accounts  for  changes  in  the  cooling  system  such  as  the  inlet 
water  temperature,  placement  temperature  of  the  concrete,  etc. 

The  command  used  to  model  the  cooling  coils  is  the  *BOUNDARY  command  in  ABAQUS.  The 
*BOUNDARY  command  is  applied  to  nodes  only.  Therefore  the  cooling  coils  need  to  be  located 
at  nodes  in  the  mesh.  When  designing  the  mesh,  it  is  important  to  consider  the  spacings  that  will 
be  attempted  for  cooling  coil  locations.  For  example,  if  the  user  would  like  to  try  4’,  6’,  8’,  and 
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10’  horizontal  spacing  of  cooling  coils,  they  may  want  to  use  a  2’  node  spacing  which  allow  for 
placement  of  cooling  coils  at  multiples  of  2’  increments. 

The  modeling  of  the  cooling  coils  employs  Y-curves  first  created  during  the  Boulder  Canyon 
Dam  Project  as  discussed  in  Chapter  1.  The  Y-curves  allow  the  ABAQUS  subroutines  to 
calculate  the  exit  temperature  of  the  cooling  coil  water  based  on  the  two  dimensionless 
parameters  given  in  Eq.  6-1  and  Eq.  6-2: 

KL  h/t 

^wPwQw  J^'2. 

(6-1)  (6-2) 


Where: 

K  =  Conductivity  of  the  concrete  to  be  cooled  in  B.t.u./ft/hr/°F 
L  =  Length  of  the  cooled  cylinder  measured  along  the  cooling  pipe  in  feet 
Ch.  =  specific  heat  of  the  cooling  water  in  B.t.u./lb/°F 
=  density  of  the  cooling  water  in  Ib/ff 
Q„  =  volume  of  water  flowing  through  the  cooling  pipe  in  ftVhr 
D  =  diameter  of  the  cooled  cylinder  in  feet 


hf  =  diffiisivity  of  the  concrete  in  ft^/hr 

t  =  total  elapsed  time  the  cooling  coil  has  been  in  operation  in  hours 


The  Y-curves  were  computed  for  a  ratio  of  b/a  of  100  where  b  is  the  radius  of  the  cooled  cylinder 
and  a  is  the  radius  of  the  cooling  pipe  as  discussed  in  Chapter  1.  Actual  cooling  pipe  spacings 
will  rarely  result  in  a  b/a  ratio  of  100.  In  order  to  take  the  difference  into  account,  charts  have 
been  tabulated  that  modify  the  values  of  h  and  D  to  account  for  this  difference.  These  modified 
values  were  shown  for  common  lift  dimensions  in  Table  2-1. 


The  Y-curves  are  shown  in  Figure  2-1.  Eq.  6-1  is  plotted  along  the  x-axis,  and  Eq.  6-2  allows  the 
user  to  locate  one  of  a  series  of  curves  plotted  on  the  Y-chart.  Using  the  intersection  point  of  the 
X- value  with  the  appropriate  curve  from  the  series,  the  Y-value  can  be  obtained.  The  Y-value  is 
related  to  the  temperatures  of  the  concrete  and  water  as  shown  in  Eq  4-3. 
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The  ABAQUS  analysis  is  broken  down  into  time  increments  within  steps.  The  model  works  by 
calculating  the  concrete  temperature  at  the  end  of  an  increment.  Using  the  Y-value  and  the  inlet 
water  temperature,  the  program  calculates  the  exit  water  temperature  for  the  increment.  The  exit 
and  inlet  water  temperatures  are  averaged  to  give  a  mean  cooling  coil  temperature  for  the 
increment.  Than  this  average  temperature  is  applied  to  the  node  where  the  cooling  coil  is  located. 

The  average  temperature  is  applied  to  the  node  at  the  desired  location  of  the  cooling  coil  using  the 
*BOUNDARY  command.  This  command  allows  the  user  to  place  a  constraint  on  a  degree  of 
freedom  in  ABAQUS.  In  the  case  of  a  cooling  coil,  a  temperature  boundary  condition  is  placed 
on  the  nodes  at  the  locations  of  the  cooling  coils.  This  boundary  temperature  is  controlled  by  a 
user  subroutine  called  DISP  which  is  called  at  the  beginning  of  each  increment  and  recalculates 
the  boundary  temperature  for  each  node  where  the  condition  is  applied.  In  this  way  the 
temperature  of  the  cooling  coil  can  be  controlled  and  kept  constant  in  the  model.  The  nodes  with 
the  temperature  boundary  condition  draw  heat  from  the  system  and  cooling  it  as  a  cooling  coil 
does  in  massive  concrete  systems. 

6.2  Implementing  the  Procedure  into  an  Abaqus  Model 

Using  the  procedure  requires  the  user  to  modify  three  things  in  the  ABAQUS  model.  These  are 
the  creation  of  node  sets,  implementing  the  *BOUNDARY  command,  and  entering  data  into  the 
user  subroutines.  The  implementation  is  exactly  the  same  for  both  2D  and  3D  models  except  for 
the  cooling  node  choices. 

6.2.1  Node  Sets  Required 


Implementation  of  this  procedure  requires  the  creation  of  several  node  sets.  All  node  sets  must 
have  the  names  specified  here  because  the  subroutines  reference  the  node  sets  by  these 
designations.  The  first  node  sets  required  are  those  which  include  the  node  numbers  of  the 
location  of  all  the  cooling  coils.  These  node  sets  have  the  designation  CCOl,  CC02,  . . .  CC12  . . . 
CCN  where  N  is  the  total  number  of  lifts.  The  numbers  01,  02,  etc.  represent  the  cooling  coils  at 
the  bottom  of  the  lift  with  the  same  number.  Therefore  CCOl  are  the  cooling  coils  between  Lift  1 
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and  the  foundation  and  CC13  are  the  cooling  coils  between  Lift  12  and  Lift  13.  Because  the 
cooling  coils  are  included  in  a  node  set  the  user  can  modify  the  number  and  location  of  the 
cooling  coils  by  changing  the  numbers  or  amount  of  the  nodes  included  in  one  of  the  CC  sets  to 
act  as  the  cooling  coils  for  that  lift. 

The  second  node  sets  required  are  necessary  to  calculate  the  concrete  temperature  at  the  end  of 
each  increment.  These  node  sets  have  the  d^ignation  AVGOl,  AVG02, ...  AVG12  ...  AVON 
where  N  is  the  total  number  of  lifts.  The  numbering  works  the  same  as  the  numbering  for  the  CC 
node  sets.  Each  ‘AVG’  node  set  should  include  two  nodes.  These  two  nodes  will  give  a 
representation  of  the  concrete  temperature.  To  locate  these  two  nodes,  the  user  should  find  the 
radius  of  the  cylinder  cooled  by  their  particular  spacing  for  coils  found  in  Table  1.  Then  33.3% 
of  the  radius  should  be  taken  as  the  distance  above  and  below  the  cooling  coil  to  find  the  average 
concrete  temperature.  The  closest  comer  node  to  this  distance  should  be  used.  Only  two  nodes, 
one  above  and  one  below  the  cooling  coils,  will  be  used  to  calculate  the  average  for  all  the  nodes 
on  that  lift.  This  is  a  valid  assumption  because  the  temperature  remains  relatively  constant 
horizontally  along  the  lift.  In  addition,  it  is  important  that  the  same  average  temperature  be  found 
for  all  the  cooling  nodes  in  one  lift  because  the  exit  temperature  of  the  water  is  dependent  on  this 
value  and  since  they  are  part  of  the  same  coil  they  must  all  have  the  same  exit  temperature.  For 
the  3D  models,  the  average  nodes  should  be  above  and  below  the  plane  of  the  cooling  coil. 

For  example,  if  the  user’s  mesh  between  Lift  01  and  Lift  02  was  as  shown  in  Figure  6-1: 

Figure  6-1:  2D  Example  Mesh 


Node  1040 


Liftl 
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In  this  case  node  set  ‘CC02’  would  include  2020,  2040,  and  2060  because  the  cooling  coils  are 
located  at  these  points  for  the  analysis.  Node  set  ‘AVG02’  would  include  1040  and  3040. 

Three  dimensional  models  are  handled  in  much  the  same  way  as  shown  in  Figure  6-2. 

Figure  6-2:  3D  Example  Mesh 


This  mesh  shows  the  location  of  the  cooling  coil  on  the  top  face  of  a  lift  with  the  lift  above  not 
shown.  If  the  user  wanted  the  cooling  coil  to  wind  in  the  location  shown,  every  bold  node  would 
need  to  be  included  in  the  node  set  'CC  corresponding  to  the  lift  number  of  the  lift  above.  Any 
cooling  coil  pattern  can  be  used  as  long  as  every  node  along  the  path  is  included  in  the  node  set 
'CC.'  The  user  also  needs  to  specify  two  average  nodes  for  node  set  'AVG.'  One  node  should  be 
chosen  above  and  below  the  plane  of  the  coils,  at  or  near  the  halfway  point  of  the  coil  length.  The 
average  nodes  chosen  should  be  corner  nodes. 

Lastly,  the  user  need  to  create  a  node  set  called  ‘ALL’  which  includes  all  the  ‘AVG’  and  ‘CC’ 
node  sets. 

The  user  should  keep  all  these  node  sets  together  in  a  section  of  the  ABAQUS  program  called 
‘Cooling  Coils.’  The  section  should  come  directly  before  the  steps,  and  the  format  is  shown 
below. 
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*************************************************  COOLING  COILS 

**************Ljp];’l.pOUNDATION 
**  COOLING  COIL  NODES 
*NSET,  NSET=CC01 
29,  33,  37 

**  AVERAGE  CONCRETE  TEMP.  NODES 
*NSET,  NSET=AVG01 
129,  5020 

»**;***********y  .TFT9-T  .TFT1 
**  COOLING  COIL  NODES 
*NSET,  NSET=CC02 
11,  15,  19 

**  AVERAGE  CONCRETE  TEMP.  NODES  ‘ 

*NSET,  NSET=AVG02 
136,  338 
** 

**************LJPP3_LjPY2 
**  COOLING  COIL  NODES 
*NSET,  NSET=CC03 
242,  246,  250 

**  AVERAGE  CONCRETE  TEMP.  NODES 
*NSET,  NSET=AVG03 

344,  538 
** 

*  4$  4: 4: ^  4: 4: 4:  *  *  *  * 

*NSET,  NSET=ALL 

CCOl,  CC02,  CC03,  AVGOl,  AVG02,  AVG03 


6.2.2  ^BOUNDARY  and  *NODE  FILE  Commands 

The  *BOUNDARY  command  needs  to  be  included  in  every  step  in  the  analysis.  The  node  sets 
for  all  the  cooling  coils  activated  in  the  step  need  to  be  included  as  shown  below. 


^BOUNDARY,  OP=NEW,  USER 
CCOl,  11,  11 
CC02,  11,  11 
CC03,  11,  11 


The  USER  command  indicates  that  the  boundary  temperature  applied  is  controlled  by  a  user 
subroutine  (DISP).  The  node  sets  are  given,  followed  by  the  starting  and  ending  degrees  of 
freedom  which  are  both  1 1  (degree  of  freedom  for  temperature). 


If  the  user  would  like  to  activate  cooling  coils  only  during  specific  steps  in  the  analysis,  this  is 
easily  accomplished  by  including  the  *BOUNDARY  command  but  not  including  specific  sets 
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within  the  command.  The  command  still  needs  to  be  included  because  the  OP=NEW  parameter 
erases  all  the  applied  boundaries  already  in  place.  This  is  shown  below. 

^BOUNDARY,  OP=NEW,  USER 

In  addition,  *NODE  FILE  must  be  included  in  every  step.  The  node  set  ‘All’  should  be  indicated 
on  the  command  line.  For  the  variable  to  be  written  to  the  results  file  the  user  should  type  ‘NT’ 
as  shown  below. 

*NODE  FILE,  NSET=ALL 
NT 

It  is  important  to  remember  not  to  activate  the  coils  for  a  particular  lift  until  the  lift  has  been  put 
in  place!  In  the  following  example,  Lift  1  is  placed  in  Step  1  and  the  cooling  coils  between  Lift  1 
and  the  Foundation  are  activated  in  this  step  (CCOl).  Lift  2  is  not  added  until  Step  4  at  which 
time  the  cooling  coils  between  Lift  1  and  Lift  2  are  activated. 


**  STEP  1,  place  lift  no.  1 
*STEP 

*HEAT  TRANSFER 
.25,2 

♦MODEL  CHANGE,  REMOVE 
START 

♦FILM,  AMPLITUDE=SUMMER,  OP=NEW 

USELS,F4„90.288 

DSELS,F2„90.288 

TLT01,F3„90.288 

♦BOUNDARY,  OP=NEW,  USER  ^ - 

CCOl,  11,  11 

♦NODE  nLE,NSET=ALL 
NT 

*RESTART,WRITE.FREQUENCY=2 

♦ENDSTEP 

**  STEP  2 
** 

♦STEP 

♦HEAT  TRANSFER 
.5,3 

♦FILM,  AMPLrrUDE=SUMMER,  OP=NEW 

USELS,F4„90.288 

DSELS,F2„90.288 

TLT01,F3„90.288 

♦BOUNDARY,  OP=NEW,  USER  ^ _ ! 

CCOl,  11,  11 

♦NODE  nLE,NSET=ALL 
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NT 

♦RESTART, WRITE, FREQUENCY=2 

♦ENDSTEP 

** 

**  STEP  3 
♦STEP 

♦HEAT  TRANSFER 

1.0,2 

♦HLM,  AMPLrnJDE=SUMMER,  OP=NEW 

USELS,F4„90.288 

DSELS,F2„90.288 

TLT01,F3„90.288 

♦BOUNDARY,  OP=NEW,  USER  ^ - 

CCOl,  11,  11 

♦NODE  nLE,NSET=ALL 
NT 

♦RESTART,WRITE,FREQUENCY=2 

♦ENDSTEP 

** 

♦♦  STEP  4,  place  lift  no.  2 
♦STEP 

♦HEAT  TRANSFER 
.25,2 

♦MODEL  CHANGE,  ADD 
LFT02 

♦FILM,  AMPLITUDE=SUMMER,  OP=NEW 

USELS,F4„90.288 

DSELS,F2„90.288 

TLT02.F3„90.288 

♦BOUNDARY,  OP=NEW,  USER  ^ - 

CCOl,  11,  11 
CC02,  11,  11 

♦NODE  nLE,NSET=ALL 
NT 

♦RESTART,  WRITE,FREQUENCY=2 
♦ENDSTEP 

ETC. 

6.2.3  Subroutine  Modification 


The  user  is  required  to  enter  certain  information  into  both  the  DISP  and  URDFE.  subroutines 
which  must  be  included  at  the  end  of  the  ABAQUS  input  file.  In  the  DISP  routine,  the  user  must 
enter  the  input  water  temperature  (I),  the  concrete  placement  temperature  (2),  and  the  Y- value  for 
the  first  increment  (3).  In  addition,  Y-values  for  every  increment  time  and  length  of  coil  in  the 
analysis  must  be  entered  (4).  As  shown  below,  the  subroutine  has  Y-values  defined  for  0.25, 0.5, 
1.0,  5.0,  7.0,  14.0,  and  15.0  day  increments  and  coil  lengths  of  480  and  600  feet.  The  user  should 
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calculate  these  Y- values  using  the  Y-curves  as  described  in  section  6.1.  In  Eq.  6-1,  the 
parameters  K,  c„,  py„  and  Qy,  are  known  and  should  not  change  during  the  analysis.  The  length  of 
the  coil  may  change  however  as  the  lifts  change  in  size.  Therefore,  different  lifts  may  have 
different  lengths  of  coil.  The  different  coil  lengths  are  associated  with  their  respective  lift 
numbers  in  the  loop  at  location  (6).  In  Eq.  6-2,  the  h,  D,  and  t  (incremental  time)  values  will 
change  with  different  cooling  coil  spacing  and  incremental  times  within  the  analysis.  The  user 
should  calculate  different  values  for  Y  based  on  these  different  increment  times.  If  a  different 
increment  time  or  length  other  than  those  listed  above  is  used,  the  user  must  add  it  at  location  (5) 
and  accommodate  it  in  the  loop  at  location  (6).  The  loop  at  location  (6)  chooses  the  correct  Y- 
value  based  on  the  increment  time  and  the  lift  number.  The  lift  number  accounts  for  the  different 
lengths  of  coil  for  different  lifts  if  necessary.  Lastly,  the  user  must  change  the  dimensions  of  the 
arrays  KCC  and  KAVG  if  necessary,  as  explained  in  number  (7).  The  DISP  routine  is  shown 
below: 


SUBROUTINE  DISP(U,  KSTEP,  KING,  TIME,  NODE,  JDOF) 

INCLUDE  ’ABA_PARAM.INC’ 

COMMON  KCC,KAVG,KTIME 
REAL*8  KAVG,KTIME 
DIMENSION  U(3),  TIME(2) 

C  A - 

C  KCC  AND  KAVG  ARE  DIMENSIONED  AS  KCC(NUMBER  OF  COIL  NODES,  NUMBER  OF 
C  LIFTS)  AND  KAVG(2,  NUMBER  OF  UFTS).  IF  THE  MODEL  REQUIRES  MORE 
C  LIFTS  OR  COILS  THAN  LISTED  CURRENTLY  SHOWN,  THE  USER  MUST  CHANGE  THE 
C  DIMENSION!  1 1  AFTER  CHANGING  THE  DIMENSIONS  FOR  KCC  AND  KAVG,  THE 
C  USER  MUST  ALSO  CHANGE  THE  ’NCOIL’  VALUE  TO  THE  MAXIMUM  NUMBER  OF 
C  COIL  NODES  PER  LIFT  AND  THE  ’NLIFT’  VALUE  TO  THE  MAXIMUM  NUMBER  OF 
C  LIFTS.  THE  DIMENSIONS  FOR  THE  ARRAYS  SHOULD  CORRESPOND  TO  THE  NCOIL 
C  AND  NLIFT  DEFINITIONS. 

C 

DIMENSION  KCC(20,50),  KAVG(2,50) 

NCOIL=20  ^ _ 

NLIFT=50 

C 

C  DEHNE  VARIABLES 
C 

C  ENTER  THE  INITIAL  COOLING  WATER  TEMPERATURE  IN  DEG.  F 

WTEMP=50.00  ◄ - 

C 

C  ENTER  THE  PLACEMENT  TEMPERATURE  OF  CONCRETE 

PT=75.00  <4 _ 

C 

C  ENTER  THE  Y-VALUE  FOR  THE  FIRST  INCREMENT 

Yn=0.22772  ^ _ 
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1 

2 

3 


on  no  on  on  on  on  on  on  on  on  on  on  on  on  on 
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C 


ENTER  THE  Y- VALUE  CORRESPONDING  TO  A  TIME  OF  0.25^DAYS  AND  280  FEET 
Y25L280=0.22772  ^ 

ENTER  THE  Y- VALUE  CORRESPONDING  TO  A  TIME  OF  0.50  DAYS  AND  280  FEET 
Y50L280=0. 19146 

ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  1.0  DAYS  AND  280  FEET 
Y1L280=0.163 

ENTER  THE  Y=VALUE  CORRESPONDING  TO  A  TIME  OF  5.0  DAYS  AND  280  FEET 
Y5L280=0.13 

ENTER  THE  Y-VALUE  CORRESPONDWG  TO  A  TIME  OF  7.0  DAYS  AND  280  FEET 
Y7L280=0.12 

ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  14.0  DAYS  AND  480  FEET 
Y14L280=0.11 

ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  15.0  DAYS  AND  480  FEET 
Y15L280=0.11 

ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  0.25  DAYS  AND  600  FEET 
Y25L600=0.24 

ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  0.50  DAYS  AND  600  FEET 
Y50L600=0.217 

ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  1.0  DAYS  AND  600  FEET 
Y1L600=0.18 

ENTER  THE  Y=VALUE  CORRESPONDING  TO  A  TIME  OF  5.0  DAYS  AND  600  FEET 
Y5L600=0.15 

ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  7.0  DAYS  AND  600  FEET 
Y7L600=0.143 

ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  14.0  DAYS  AND  600  FEET 
Y14L600=0.12 

ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  15.0  DAYS  AND  600  FEET 
Y15L600=0.12  ^ - 
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ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF 
Y  = 


_  DAYS  AND _ FEET 

◄ -  5 

Added  if  necessary  by  user 


IF(KSTEP.EQ.l  .AND.  KINC.EQ.1)THEN 
TINC=TFI 
ELSE 

TINC=TIME(2)-KTIME 
END  IF 

DO  50  N=l,NCOIL 
DO  40  M=1,NLIFT 


IF  (NODE.EQ.KCC(N,M))  THEN 
LIFT=M 
END  IF 

40  CONTINUE 
50  CONTINUE 
C 

IF  (TINC.LE.0.25  .AND.  LIFT.EQ.  I )  THEN 

Y=Y25L600  ^ - 

ELSE  IF  (TINC.GT.0.25  .AND.  TINC.LE.0.5  .AND.  LIFT.EQ.  1)  THEN 
Y=Y50L600 

ELSE  IF  (TINC.GT.0.5  .AND.  UNCLE.  1.0  .AND.  LIFT.EQ  1)  THEN 
Y=Y1L600 

ELSE  IF  (TINC.GT.  1.0  .AND.  TINC.LE.5.0  .AND.  LIFT.EQ.  1)  THEN 
Y=Y5L600 

ELSE  IF  (TINC.GT.5.0  .AND.  TINC.LE.7.0  .AND.  LIFT.EQ  1)  THEN 
Y=Y7L600 

ELSE  IF  (TINC.GT.7.0  .AND.  TINC.LE.  14.0  .AND.  LIFT.EQ.  1 )  THEN 
Y=Y14L600 

ELSE  IF(TINC.EQ.15.0  .AND.  LIFT.EQ.1)  THEN 
Y=Y15L600 

ELSE  IF(TINC.LE.0.25  .AND.  LIFT.EQ.2)  THEN 
Y=Y25L280 

ELSE  IF  (TINC.GT.0.25  .AND.  TINC.LE.0.5  .AND.  LIFT.EQ.2)  THEN 
Y=Y50L280 

ELSE  IF  (TINC.GT.0.5  .AND.  TINC.LE.  1.0  .AND.  LIFT.EQ.2)  THEN 
Y=YIL280 

ELSE  IF  (TINC.GT.1.0  .AND.  TINC.LE.5.0  .AND.  LIFT.EQ.2)  THEN 
Y=Y5L280 

ELSE  IF  (TINC.GT.5.0  .AND.  TINC.LE.7.0  .AND.  LIFT.EQ  2)  THEN 
Y=Y7L280 

ELSE  IF  (TINC.GT.7.0  .AND.  TINC.LE.  14.0  .AND.  LIFT  EQ  2)  THEN 
Y=Y14L280 

ELSE  IF  (TINC.EQ.15.0  .AND.  LIFT.EQ.2)  THEN 
Y=Y15L280 

ELSE 

Y=Y15 

END  IF 


CONCT=(KAVG(  1  ,LIFT)+KAVG(2,LIFT))/2 
C 

IF  (KINC.EQ.l  .AND.  KSTEP.EQ.1)  THEN 

U(  1  )=((YFI*(PT-WTEMP)+WTEMP)+WTEMP)/2 
ELSE 

U(  1  )=((Y*(CONCT-WTEMP)+WTEMP)+WTEMP)/2 
END  IF 
C 

RETURN 

END 


This  loop  is 
modifled  as 
necessary  by  the 
user.  Create  a 
position  for  the 
time  increment, 
lift  number,  and 
the  Y- value  for 
that  increment 
time  and  length 
of  coil  for  that 
narticular  lift. 
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There  is  only  one  parameter  in  addition  to  the  array  dimensions  discussed  above  that  needs  to  be 
modified  in  the  subroutine  URDFIL.  This  is  the  file  name.  The  file  name  that  needs  to  be 
entered  is  the  name  of  the  input  file  for  the  particular  model  being  run.  It  is  entered,  without  the 
ending,  in  single  quotation  marks.  For  the  example  shown  below,  the  file  name  was  modelr.inp 
(8).  Lastly,  the  user  must  change  the  dimensions  of  the  arrays  KCC  and  KAVG  if  necessary,  as 
explained  in  number  (7)  above  for  the  DISP  routine.  Only  the  beginning  of  the  URDFIL 
subroutine  is  shown  below  because  this  is  the  only  portion  of  the  routine  which  requires 
modification  by  the  user. 


SUBROUTINE  URDnL(LSTOP,LOVRWRT,KSTEP,KINC) 

INCLUDE  ’ABA_PARAM.INC’ 

C  ◄ - 1 

C  KCC  AND  KAVG  ARE  DIMENSIONED  AS  KCC(NUMBER  OF  COIL  NODES,  NUMBER  OF 
C  LIFTS)  AND  KAVG(2,  NUMBER  OF  LIFTS).  IF  THE  MODEL  REQUIRES  MORE 
C  LIFTS  OR  COILS  THAN  LISTED  CURRENTLY  SHOWN,  THE  USER  MUST  CHANGE  THE 
C  DIMENSION! ! !  AFTER  CHANGING  THE  DIMENSIONS  FOR  KCC  AND  KAVG,  THE 
C  USER  MUST  ALSO  CHANGE  THE  ’NCOIL’  VALUE  TO  THE  MAXIMUM  NUMBER  OF  _ 

C  COIL  NODES  PER  LIFT  AND  THE  MJFT’  VALUE  TO  THE  MAXIMUM  NUMBER  OF 
C  LIFTS.  THE  DIMENSIONS  FOR  THE  ARRAYS  SHOULD  CORRESPOND  TO  THE  NCOIL 
C  AND  NLIFT  DEFINITIONS. 

C 

DIMENSION  KCC(20,50),IAVG(2,50),  KAVG(2,50) 

DIMENSION  ARRAY(513),JRRAY(NPRECD,513) 

EQUIVALENCE  (ARRAY(1),JRRAY(1,1)) 

COMMON  KCC,KAVG,KTIME 
CHARACTER*80  FNAME 
REAL*8  KAVG,KTIME 
NCOIL=20 

NLIFT=50  ◄ - - - - - 

ICOUNT1=0 

ICOUNT2=0 

IVAR=0 

C 

C  ENTER  THE  NAME  OF  THE  INPUT  HLE  IN  SINGLE  QUOTES  IWTHOUT  THE 
C  EXTENTION 

FNAME=’modelr’  - -  8 

C 

DO  100  Kl=l,99999 
C 

DO  101=1,513 
ARRAYa)=0.0 
10  CONTINUE 
C 


For  the  user  who  wants  a  more  in  depth  explanation  of  the  subroutines,  their  explanation  is 
included  in  the  appendix. 
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Chapter  7 
Conclusions 


This  chapter  will  summarize  the  advantages  and  disadvantages  of  the  boundary  temperature 
modeling  procedure.  It  will  also  discuss  possible  future  work. 

7.1  Model  Advantages  and  Disadvantages 


The  boundary  node  method  for  modeling  the  cooling  coils  has  several  advantages,  and  also 
disadvantages  as  outlined  below. 

7.1.1  Advantages 


A  major  advantage  of  using  the  boundary  method  of  modeling  the  cooling  coils  is  that  it  is  easily 
implemented  into  an  existing  model.  No  change  of  the  mesh  is  required,  and  no  special  mesh, 
element  size,  or  node  placement  is  required  to  use  the  procedure.  The  actual  mechanics  of 
implementing  the  procedure  into  a  model  are  outlined  in  Chapter  6.  However,  the 
implementation  is  not  rigorous,  and  requires  no  subroutine  modification  by  the  user  beyond 
entering  certain  parameters  such  as  Y-values,  placement  concrete  temperature,  and  input  water 
temperature.  Also,  implementation  is  essentially  the  same  in  the  two-dimensional  and  three- 
dimensional  analysis,  reducing  confusion  for  the  user. 
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The  two-dimensional  model  would  be  most  likely  used  to  try  a  variety  of  cooling  coil  spacings 
and  designs.  It  is  easy  to  modify  the  cooling  coil  spacings  as  long  as  the  user  has  nodes  placed 
where  they  would  like  to  place  coils.  In  addition,  the  three  dimensional  model  can  accommodate 
any  winding  pattern  for  the  cooling  coils  layout. 

Another  important  point  to  emphasize  is  that  the  model  does  not  require  that  the  cooling  coils  be 
functional  for  the  duration  of  the  analysis.  The  coils  can  be  turned  on  for  any  steps  in  the 
analysis,  then  turned  off  when  the  user  would  like  to  stop  cooling.  They  can  also  be  turned  back 
on  later  in  the  analysis  if  desired  by  the  user. 

7.1.2  Disadvantages 

The  most  important  disadvantage  to  the  boundary  node  method  of  modeling  the  cooling  coils  is 
the  reliance  on  the  Y-curves.  The  Y-curves  were  developed  based  on  several  assumptions, 
including  an  insulated  cylinder,  a  straight  coil,  and  non-incremental  analysis.  However,  the 
curves  have  been  used  extensively  by  the  Army  Corp  of  Engineers  since  the  construction  of  the 
Hoover  Dam,  with  good  results.  In  addition,  as  has  been  established,  the  error  due  to  an 
incremental  analysis  is  minor.  However,  a  method  of  analysis  which  does  not  require  the  use  of 
the  Y-curves  would  be  desirable. 

The  other  disadvantage  of  using  boundary  nodes  to  model  the  coils  is  that  the  coils  must  be 
placed  at  node  locations.  Due  to  this  constraint,  the  user  should  consider  possible  spacing  options 
before  creating  their  mesh  in  order  to  choose  a  node  spacing  which  allows  for  analysis  of  several 
coil  spacing  options. 

7.2  Future  Work 

The  development  of  a  modeling  procedure  that  could  be  used  to  verify  the  Y-curves  is  still 
sought-after.  It  does  not  necessarily  have  to  be  implemented  for  design,  only  to  verify  the  curves. 
The  convection/diffusion  modeling  methods  available  to  model  forced  convection  through  a  mesh 
in  ABAQUS,  while  proven  not  to  be  adequate  for  design,  may  still  have  the  potential  to  be  used 
as  a  verification  of  the  Y-curves. 
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Appendix  A 

Subroutine  URDFIL 


SUBROUTINE  URDnL(LSTOP,LOVRWRT,KSTEP,KINC) 

INCLUDE  ’ABA.PARAM.INC’ 

C 

C  KCC  AND  KAVG  ARE  DIMENSIONED  AS  KCC(NUMBER  OF  COIL  NODES,  NUMBER  OF 
C  LIFTS)  AND  KAVG(2,  NUMBER  OF  LIFTS).  IF  THE  MODEL  REQUIRES  MORE 
C  LIFTS  OR  COILS  THAN  LISTED  CURRENTLY  SHOWN,  THE  USER  MUST  CHANGE  THE 
C  DIMENSION!!!  AFTER  CHANGING  THE  DIMENSIONS  FOR  KCC  AND  KAVG  THE 
C  USER  MUST  ALSO  CHANGE  THE  NCOIL’  VALUE  TO  THE  MAXIMUM  NUMBER  OF 
C  COIL  NODES  PER  LIFT  AND  THE  ’NLIFT’  VALUE  TO  THE  MAXIMUM  NUMBER  OF 
C  LIFTS.  THE 

C  DIMENSIONS  FOR  THE  ARRAYS  SHOULD  CORRESPOND  TO  THE  NCOIL  AND  NLIFT 
C  DEFINITIONS. 

C 

DIMENSION  KCC(20,50),IAVG(2,50),  KAVG(2,50) 

KCC  is  the  array  which  stores  the  node  numbers  of  the  cooling  coil  nodes.  It  assumes  no  more 
than  20  cooling  coils  per  lift  and  50  lifts  unless  modified  as  stated  above.  lAVG  is  an  array 
which  stores  the  node  numbers  of  the  average  concrete  temperature  nodes.  There  are  only 
two  average  nodes  per  lift  and  it  again  assumes  50  lifts  unless  modified.  KAVG  is  an  array 
which  stores  the  temperature  of  the  nodes  in  lAVG.  NCOIL  and  NLIFT  correspond  to  the 
maximum  number  of  coils  and  lifts  as  defined  in  the  dimension  statements  for  KCC  KAVG 
and  lAVG. 

DIMENSION  ARRAY(513),JRRAY(NPRECD,513) 

EQUIVALENCE  (ARRAY(  I  ).JRRAY(  1,1)) 

COMMON  KCC,KAVG,KTIME 

These  three  variables  need  to  be  passed  to  the  subroutine  DISP  and  are  called  out  as  common 
blocks. 

CHARACTER*80  FNAME 

REAL*8  KAVG,KTIME 

NCOIL=20 

NLIFT=50 

ICOUNT1=0 

ICOUNT2=0 
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IVAR=0 

C 

C  EOTER  THE  NAME  OF  THE  INPUT  FILE  IN  SINGLE  QUOTES  WITHOUT  THE 

C  EXTENTION 

FNAME=Ynodelr’  . 

File  name  of  the  *.fil  (results  file)  which  Abaqus  reads  to  get  the  node  temperature  during  the 

analysis. 

C 

DO  100  Kl=l,99999 

C 

DO  101=1,513 

ARRAY(I)=0.0 
10  CONTINUE 

This  loop  initializes  the  values  of  Array  to  0. 

C 

CALL  DBFILE(0,ARRAY,JRCD) 

IF(JRCD.NE.0)GO  TO  1 10 

KEY=JRRAY(1,2)  .  . 

Every  time  this  call  to  DBFBLE  is  made  Abaqus  calls  a  new  record.  The  following  code  is 
completed  for  that  particular  record  which  may  or  may  not  have  the  information  needed. 
Then  the  next  record  is  called  and  the  process  is  repeated  until  all  the  records  have  been  read. 
When  this  happens  JRCD  is  set  to  something  other  than  zero  and  the  routine  ends.  KEY  is 
the  number  designation  of  every  record. 

C 

IF  (KEY.EQ.  1931)  THEN 

KEY  1931  contains  the  information  for  node  sets  Including  the  node  set  name  and  the 
numbers  of  all  the  nodes  included  in  the  set. 

C 

IF  (ARRAY(3).EQ.CC01’)  THEN 
ICOUNTl=l 

ELSE  IF  (ARRAY(3).EQ.CC02’)  THEN 
ICOUNTl=2 

ELSE  IF  (ARRAY(3).EQ. Ceos’)  THEN 
ICOUNTl=3 

ELSE  IF  (ARRAY(3).EQ.CC04’)  THEN 
ICOUNTl=4 

ELSE  IF  (ARRAY(3).EQ.CC05’)  THEN 
ICOUNTl=5 

ELSE  IF  (ARRAY(3).EQ.CC06’)  THEN 
ICOUNTl=6 

ELSE  IF  (ARRAY(3).EQ.CC07’)  THEN 
ICOUNTl=7 

ELSE  IF  (ARRAY(3).EQ. CCDS’)  THEN 
ICOUNTl=8 

ELSE  IF  (ARRAY(3).EQ.CC09’)  THEN 
ICOUNTl=9 

ELSE  IF  (ARRAY(3).EQ.CC10’)  THEN 
ICOUNT1=10 

ELSE  IF  (ARRAY(3).EQ.CC1  H  THEN 
ICOUNTl=ll 

ELSE  IF  (ARRAY(3).EQ.CC12’)  THEN 
ICOUNTl=12 

ELSE  IF  (ARRAY(3).EQ. eels’)  THEN 
ICOUNTl=13 


o  n 
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ELSE  IF  (ARRAY(3).EQ.X:C14’)  THEN 
IC0UNT1=14 

ELSE  IF  (ARRAY(3).EQ. CGIS’)  THEN 
IC0UNT1=15 

ELSE  IF  (ARRAY(3).EQ.CC16’)  THEN 
IC0UNT1=16 

ELSE  IF  (ARRAY(3).EQ.CC17’)  THEN 
ICOUNTl=17 

ELSE  IF  (ARRAY(3).EQ. CGIS’)  THEN 
IG0UNT1=18 

ELSE  IF  (ARRAY(3).EQ.CG19’)  THEN 
IGOUNTl=19 

ELSE  IF  (ARRAY(3).EQ.CG20’)  THEN 
IGOUNT1=20 

ELSE  IF  (ARRAY(3).EQ.CG21’)  THEN 
IG0UNT1=21 

ELSE  IF  (ARRAY(3).EQ.CG22’)  THEN 
IG0UNT1=22 

ELSE  IF  (ARRAY(3).EQ.CG23’)  THEN 
IG0UNT1=23 

ELSE  IF  (ARRAY(3).EQ.CG24’)  THEN 
IG0UNT1=24 

ELSE  IF  (ARRAY(3).EQ.CG25’)  THEN 
IG0UNT1=25 

ELSE  IF  (ARRAY(3).EQ.CG26’)  THEN 
IG0UNT1=26 

ELSE  IF  (ARRAY(3).EQ.CG27’)  THEN 
IG0UNT1=27 

ELSE  IF  (ARRAY(3).EQ.CG28’)  THEN 
IG0UNT1=:28 

ELSE 

IGOUNT1=0 

END  IF 

If  the  KEY=1931,  Array(3)  contains  the  nodeset  name.  If  Array(3)  matches  the  GG*  name 
shown  above,  the  variable  IGOUNTl  is  set  to  the  lift  number. 

C 

IF  (IGOUNTl  .GT.O)  THEN 
DO  20  N=l,NGOIL 

KGG(N,IGOUNT  1  )=JRRA  Y(  1 ,3+N) 

20  GONTINUE 

END  IF 

If  IGOUNTl  is  equal  to  a  lift  number  (greater  than  0)  then  this  loop  fills  the  array  KGG 
with  the  node  numbers  included  in  the  node  set  corresponding  to  that  lift  These  numbers 
are  in  JRRAY(l,3...End  of  set). 


IF  (ARRAYISI.EQ.’AVGOn  THEN 
IGOUNT2=l 

ELSE  IF  (ARRAY(3).EQ.’AVG02’)  THEN 
IGOUNT2=2 

ELSE  IF  (ARRAY(3).EQ.’AVG03’)  THEN 
IGOUNT2=3 

ELSE  IF  (ARRAY(3).EQ.’AVG04’)  THEN 
IGOUNT2=4 


ELSE  IF  (ARRAY(3).EQ.’AVG05’)  THEN 
ICOUNT2=5 

ELSE  IF  (ARRAY(3).EQ.’AVG06’)  THEN 
IC01INT2=6 

ELSE  IF  (ARRAY(3).EQ.’AVG07’)  THEN 
ICOUNT2=7 

ELSE  IF  (ARRAY(3).EQ.’AVG08’)  THEN 
ICOUNT2=8 

ELSE  IF  (ARRAY(3).EQ.’AVG09’)  THEN 
ICOUNT2=9 

ELSE  IF  (ARRAY(3).EQ.’AVG10’)  THEN 
ICOUNT2=10 

ELSE  IF  (ARRAY(3).EQ.’AVG1 1 0  THEN 
ICOUNT2=ll 

ELSE  IF  (ARRAY(3).EQ.’AVG12’)  THEN 
ICOUNT2=12 

ELSE  IF  (ARRAY(3).EQ.’AVG13’)  THEN 
ICOUNT2=13 

ELSE  IF  (ARRAY(3).EQ.’AVG14’)  THEN 
ICOUNT2=14 

ELSE  IF  (ARRAY(3).EQ.’AVG15’)  THEN 
ICOUNT2=15 

ELSE  IF  (ARRAY(3).EQ.’AVG16’)  THEN 
ICOUNT2=16 

ELSE  IF  (ARRAY(3).EQ.’AVG17’)  THEN 
ICOUNT2=17 

ELSE  IF  (ARRAY(3).EQ.’AVG18’)  THEN 
IC0UNT2=18 

ELSE  IF  (ARRAY(3).EQ.’AVG19’)  THEN 
ICOUNT2=19 

ELSE  IF  (ARRAY(3).EQ.’AVG20’)  THEN 
ICOUNT2=20 

ELSE  IF  (ARRAY(3).EQ.’AVG21 0  THEN 
ICOUNT2=21 

ELSE  IF  (ARRAY(3).EQ.’AVG22’)  THEN 
ICOUNT2=22 

ELSE  IF  (ARRAY(3).EQ.’AVG23’)  THEN 
ICOUNT2=23 

ELSE  IF  (ARRAY(3).EQ.’AVG24’)  THEN 
ICOUNT2=24 

ELSE  IF  (ARRAY(3).EQ.’AVG25’)  THEN 
ICOUNT2=25 

ELSE  IF  (ARRAY(3).EQ.’AVG26’)  THEN 
ICOUNT2=26 

ELSE  IF  (ARRAY(3).EQ.’AVG27’)  THEN 
ICOUNT2=27 

ELSE  IF  (ARRAY(3).EQ.’AVG28’)  THEN 
ICOUNT2=28 

ELSE 

ICOUNT2=0 

END  IF 

IF  (ICOUNT2.GT.O)  THEN 
DO  30  N=l,2 
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IAVG(N,IC0UNT2)=JRRAY(1,3+N) 

30  CONTINUE 

END  IF 
END  IF 

This  section  stores  the  node  numbers  included  in  the  KAVG*  node  sets  in  the  same  way  as 
those  above  for  the  CC*  node  sets. 

C 

IF(KEY.EQ.2000)THEN 

KEY  2000  stores  information  including  the  step  number,  increment  number,  and  total 
time. 

KISTEP=JRRAY(1,8) 

KSINC=JRRAY(1,9) 

It  is  necessary  to  know  the  total  time  when  this  routine  is  called  at  the  end  of  the 
increment  To  find  the  total  time  we  needed  to  find  the  current  step  and  increment  in  the 
*.fil  file.  To  do  this  we  needed  to  compare  the  current  step  and  increment  (read  in  as 
KSTEP  and  KING  in  the  subroutine  call)  and  compare  these  with  the  KISTEP  and  KSINC 
values  defined  above  which  give  the  step  and  increment  in  the  *.fU  file.  The  Abaqus  utility 
routine  POSFTL  could  not  be  used  to  position  us  at  the  current  step  and  increment  in  the 
*.fil  file  because  it  positioned  us  after  the  keys  1931  and  2000  which  we  needed. 

IF  (KINC.EQ.KSINC  .AND.  KSTEP.EQ.KISTEP)  TOEN 
IVAR=1 

KTIME=ARRAY(3) 

ELSE 

IVAR=0 
END  IF 
END  IF 

KTIME  is  the  total  time  at  the  end  of  the  current  increment.  IVAR  is  a  dummy  variable  to 
indicate  whether  or  not  we  are  at  the  current  step  and  increment  in  the  *.fil  file. 

C 

IF(IVAR.EQ.l  .AND.  KEY.EQ.201)THEN 

If  IVAR=1  than  we  are  positioned  at  the  current  step  and  increment  in  the  *.fil  file.  KEY 
201  contains  the  node  number  in  position  JRRAY(1,3)  and  the  node  temperature  in 
ARRAY(1,4). 

DO  60N=1,2 

DO  50M=1,NLIFT 

IF  (JRRAY(  1 ,3).EQ.IAVG(N,M))  THEN 
KAVG(N,M)=ARRAY(4)  ’ 

END  IF 

This  loop  fills  array  KAVG  with  the  node  temperatures  of  the  nodes  in  array  lAVG. 

50  CONTINUE 

60  CONTINUE 
END  IF 
C 

100  CONTINUE 
no  CONTINUE 
RETURN 
END 


noon 
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Appendix  B 


Subroutine  DISP 


SUBROUTINE  DISP(U,  KSTEP,  KING,  TIME,  NODE,  JDOF) 

INCLUDE  ’ABA_PARAM.INC’ 

COMMON  KCC,KAVG,KTIME 

These  arrays  are  read  in  from  the  URDFIL  subroutine.  The  URDFIL  subroutine  is  called  at 
the  end  of  the  increment.  DISP  is  called  at  the  beginning  of  the  increment  Therefore  KAVG, 
and  KTIME  are  the  average  node  temperatures  and  the  total  time  for  the  previous  increment 
REAL*8  KAVG,KTIME 
DIMENSION  U(3),  TIME(2) 

C 

C  KCC  AND  KAVG  ARE  DIMENSIONED  AS  KCC(NUMBER  OF  COIL  NODES,  NUMBER  OF 
C  LIFTS)  AND  KAVG(2,  NUMBER  OF  LIFTS).  IF  THE  MODEL  REQUIRES  MORE 
C  LIFTS  OR  COILS  THAN  USTED  CURRENTLY  SHOWN,  THE  USER  MUST  CHANGE  THE 
C  DIMENSION!!!  AFTER  CHANGING  THE  DIMENSIONS  FOR  KCC  AND  KAVG,  THE 
C  USER  MUST  ALSO  CHANGE  THE  ’NCOIL’  VALUE  TO  THE  MAXIMUM  NUMBER  OF 
C  COILS  PER  LIFT  AND  THE  ’NLIFT’  VALUE  TO  THE  MAXIMUM  NUMBER  OF  LIFTS. 

C  THE 

C  DIMENSIONS  FOR  THE  ARRAYS  SHOULD  CORRESPOND  TO  THE  NCOIL  AND  NLIFT 
C  DEHNIHONS. 

C 

DIMENSION  KCC(20,50),  KAVG(2,50) 

NCOIL=20 

NLIFT=50 

KCC  is  the  array  which  stores  the  node  numbers  of  the  cooling  coil  nodes.  It  assumes  no  more 
than  20  cooling  coils  per  lift  and  50  lifts  unless  modified  as  stated  above.  There  are  only  two 
average  nodes  per  lift  and  it  again  assumes  50  lifts  unless  modified.  KAVG  is  an  array  which 
stores  the  temperature  of  the  nodes  in  lAVG.  NCOIL  and  NLIFT  correspond  to  the 
maximum  number  of  coils  and  lifts  as  defined  in  the  dimension  statements  for  KCC,  KAVG, 
and  lAVG. 

DEFINE  VARIABLES 

ENTER  THE  INITIAL  COOLING  WATER  TEMPERATURE  IN  DEG.  F 
WTEMP=50.00 
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C 

C  ENTER  THE  PLACEMENT  TEMPERATURE  OF  CONCRETE 
PT=75.00 
C 

C  ENTER  THE  Y-VALUE  FOR  THE  FIRST  INCREMENT 
Yn=0.25 
C 

C  ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  0.25  DAYS  AND  480 
C  FEET 

Y25L480=0.24 

C 

C  ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  0.50  DAYS  AND  480 
C  FEET 

Y50L480=0.217 

C 

C  ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  1 .0  DAYS  AND  480 
C  FEET 

Y1L480=0.18 

C 

C  ENTER  THE  Y=VALUE  CORRESPONDING  TO  A  TIME  OF  5.0  DAYS  AND  480 
C  FEET 

Y5L480=0.15 

C 

C  ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  7.0  DAYS  AND  480 
C  FEET 

Y7L480=0.143 

C 

C  ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  1 4.0  DA YS  AND  480 
C  FEET 

Y14L480=0.12 

C 

C  ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  1 5.0  DAYS  AND  480 
C  FEET 

Y15L480=0.12 

C 

C  ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  0.25  DAYS  AND  600 
C  FEET 

Y25L600=0.24 

C 

C  ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  0.50  DAYS  AND  600 
C  FEET 

Y50L600=0.217 

C 

C  ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  1 .0  DAYS  AND  600 
C  FEET 

Y1L600=0.18 

C 

C  ENTER  THE  Y=VALUE  CORRESPONDING  TO  A  TIME  OF  5.0  DAYS  AND  600 
C  FEET 

Y5L600=0.15 

C 

C  ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  7.0  DAYS  AND  600 
C  FEET 

Y7L600=0.143 
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C 

C  ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  14.0  DAYS  AND  600 
C  FEET 

Y14L600=0.12 

C 

C  ENTER  THE  Y-VALUE  CORRESPONDING  TO  A  TIME  OF  15.0  DAYS  AND  600 
C  FEET 

Y15L600=0.12 

C 

User  deflned  variables. 

C 

IF  (KSTEP.EQ.l  .AND.  KINC.EQ.l)  THEN 
TINC=TFI 
ELSE 

TINC=TIME(2)-KTIME 
END  IF 

If  the  current  step  and  increment  are  both  one  (KSTEP  and  KINC=1),  it  is  the  first  increment 
of  the  analysis.  If  it  is  a  later  step,  TIME(2)  is  the  total  time  at  the  end  of  the  current 
increment.  By  subtracting  the  total  time  from  the  previous  increment  (KTIME),  the  increment 
time  (TINC)  is  determined. 

C 

DO  50  N=l,NCOIL 
DO  40M=1,NLIFT 

IF  (NODE.EQ.KCC(N,M))  THEN 
LIFT=M 
END  IF 

40  CONTINUE 
50  CONTINUE 

NODE  is  the  boundary  node  which  the  routine  is  being  called  for  (the  cooling  coil  that  the 
routine  is  being  called  for).  The  nodes  position  in  the  KCC  array  determines  which  lift  it 
corresponds  to  and  this  number  is  LIFT. 

C 

IF  (TINC.LE.0.25  .AND.  LIFT.EQ.l)  THEN 
Y=Y25L600 

ELSE  IF  (TINC.GT.0.25  .AND.  TINC.LE.0.5  .AND.  LIFT.EQ.l)  THEN 
Y=Y50L600 

ELSE  IF  (nNC.GT.O.S  .AND.  TINC.LE.1.0  .AND.  LIFT.EQ.l)  THEN 
Y=YIL600 

ELSE  IF  (TINC.GT.1.0  .AND.  TINC.LE.5.0  .AND.  LIFT.EQ.l)  THEN 
Y=Y5L600 

ELSE  IF  (TINC.GT.5.0  .AND.  TINC.LE.7.0  .AND.  LIFT.EQ.l)  THEN 
Y=Y7L600 

ELSE  IF  (TINC.GT.7.0  .AND.  TINC.LE.14.0  .AND.  LIFT.EQ.l)  THEN 
Y=Y14L600 

ELSE  IFaiNC.EQ.15.0  .AND.  LIFT.EQ.l)  THEN 
Y=Y15L600 

ELSE  IF(TINC.LE.0.25  .AND.  LIFT.EQ.2)  THEN 
Y=Y25L480 

ELSE  IF  (TINC.GT.0.25  .AND.  TINC.LE.0.5  .AND.  LIFT.EQ.2)  THEN 
Y=Y50L480 

ELSE  IF  (TINC.GT.0.5  .AND.  TINC.LE.L0  .AND.  LIFT.EQ.2)  THEN 
Y=Y1L480 

ELSE  IF  (TINC.GT.1.0  .AND.  TINC.LE.5.0  .AND.  LIFT.EQ.2)  THEN 
Y=Y5L480 


ELSE  IF  (TINC.GT.5.0  .AND.  TINC.LE.7.0  .AND.  LIFT.EQ.2)  THEN 
Y=Y7L480 

ELSE  IF  (TINC.GT.7.0  .AND.  TINC.LE.14.0  .AND.  LIFT.EQ.2)  THEN 
Y=Y14L480 

ELSE  IF  (TINC.EQ.  15.0  .AND.  LIFT.EQ.2)  THEN 
Y=Y15L480 
ELSE 
Y=Y15 
END  IF 

This  loop  determines  what  Y-values  should  be  used  based  on  the  current  increment  time  and 
the  lift  number.  The  reason  this  lift  number  is  important  is  because  different  lifts  have 
different  lengths  of  coil  which  affects  the  Y-value.  The  last  else  statement  is  used  as  a  default 
in  case  a  time  increment  is  used  during  the  analysis  that  does  not  have  a  Y-value.  This  should 
never  happen  unless  the  user  makes  a  mistake  in  implementing  the  procedure. 

CONCT=(KAVG(l,LIFT)+KAVG(2,LIFT))/2 

The  average  concrete  temperature  for  the  boundary  node  considered  is  CONCT. 

IF  (KINC.EQ.1  .AND.  KSTEP.EQ.1)  THEN 

U(  1  )=((Yn*(PT-WTEMP)+WTEMP)+WIEMP)/2 
ELSE 

U(  1  )=((Y*(CONCT-WTEMP)+WTEMP)+WTEMP)/2 
END  IF 

U(l)  is  the  boundary  temperature  applied  to  the  boundary  node.  If  the  analysis  is  at  the 
beginning  (increment  and  step  1),  then  the  boundary  temperature  is  calculated  using  the 
placement  temperature  of  the  concrete.  If  it  is  later  in  the  analysis,  the  current  concrete 
temperature  (CONCT)  is  used  instead. 

RETURN 

END 
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